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INTRODUCTION 


During  1982-P3  work  progressed  on  four  separate  fronts 
corresponding  to  both  experimental  and  computational  aspects  of 
turbulent  flow  in  1P0°  bends  for  both  square  and  round  tubes. 

The  year  saw  the  completion  of  work  on  the  square  sectioned'' 
bend  and  Mr.  R.W.  Johnson's  PhD  thesis  documenting  in  great  detail 
the  research  accomplished  will  be  circulated  within  two  months. 
Here,  therefore,  in  Section  2  only  a  short  summary  will  be  given. 
Appendix  1  includes  a  paper  by  Chang,  Humphrey,  Johnson  and 
Launder  [l]  reporting  the  outcome  of  the  computational  work 
reported  at  the  4th  Turbulent  Shear  Flows  Symposium  in  Karlsruhe. 


The  computational  work  on  flow  in  circular  sectioned  ducts 
has  also  reached  the  stage  of  publication  [2,3,4]  and  copies  of 
these  papers  are  also  appended  to  supplement  the  summary  of 
computational  and  experimental  research  given  in  Section  3. 
Finally,  Section  4  outlines  the  work  now  underway  to  bring  the 
project  to  completion.  / 


SQUARE  SECTIONED  BEND 
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>aratus  and  Instrumentation 


2.1  Experimental  Program  ^  A 

X 

2.1.1  Apparatus  and  Instrumentation  The  design  and  construction 
of  the  apparatus  was  described  in  some  detail  in  last  year's 
Annual  Report  464^*  The  radius  of  the  bend  was  3.35  times  the  duct 
hydraulic  diameter  (D^> ,  the  width  of  the  duct  walls  being  RP.9 
mm.  Up  to  70  hydraulic  diameters  of  flow  development  were 
available  prior  to  the  bend  but  the  main  emphasis  of  the 
experiments  was  directed  at  the  case  of  31  Dg  which  corresponded 
with  the  conditions  at  which  the  UC  Berkeley  flow  field  data  had 
been  taken.^(In  fact,  no  major  differences  in  tbe^flow  or  thermal 
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behaviour  in  the  bend  i'tself  arose  from  the  variation  in  the 
length  of  development) .  ^Downstream  of  the  bend  40  hydraulic 
diameters  of  instrumented  straight  ducting  allowed  the  effects  of 
decaying  secondary  flow  to  be  studied. 


Heating  was  provided  by  electrically-heated  Intrex 
sheeting,  a  gold  film  deposited  uniformly  on  a  plastic  substrate. 
The  sheets  were  cut  to  size  and  affixed  to  the  inside  of  the  duct 
with  the  gold  film  on  the  surface  exposed  to  the  airstream.- 
Variations  in  resistivity  of  the  Intrex  were  found  to  be  within  - 
8%  of  the  mean  provided  material  from  near  the  edges  of  the  roll 
was  discarded.  In  use  a  certain  amount  of  ageing  was  evident  and 
when  local  non-uniformities  became  serious  the  sheets  were  removed 
and  replaced  by  new  material. 


Thermocouple  measurements  at  up  to  9  stations  recorded  wall 
temperatures  at  12  positions  around  one  half  of  the  perimeter 
(the  flow  being  nominally  symmetric  about  the  surface  passing 
through  the  centre  of  the  duct  and  lying  in  the  plane  of  the 
bend) .  Temperatures  were  also  measured  at  3  positions  on  the 
opposite  side  of  the  duct  to  provide  symmetry  checks.  The  maximum 
wall  temperature  was  limited  to  30°C  above  ambient  to  limit  the 
rate  of  deterioration  of  the  Intrex  sheeting  and  to  keep  the 
variation  in  air  properties  to  unimportant  levels.  A  rake 
carrying  13  chromel-alumel  thermocouples  allowed  measurements  of 
the  temperature  field  over  the  duct  cross  section  at  the  same 
positions  as  the  wall  temperatures  were  recorded. 

The  ducting  was  encased  in  expanded  polystyrene  insulant  to 
a  thickness  of  40  mm. 


2.1.2  Test  Program  Measurements  of  wall  and  interior 
temperatures  and  heat  fluxes  were  made  at  a  Reynolds  number  of 
58700  corresponding  to  that  of  the  Berkeley  experiment.  Moreover, 
although  Professor  Humphrey's  group  were  responsible  for 


documenting  the  flow  field  [fi] ,  some  limited  velocity-field 
measurements  were  judged  necessary  at  UMIST  to  establish  that  our 
flow  conditions  were  indeed  sensibly  the  same  as  at  Berkeley  (this 
verification  acquired  further  importance  as  it  was  found  that  the 
numerical  computations  exhibited  a  strikingly  different  flow 
pattern  midway  around  the  bend  from  the  measurements) .  Velocity 
profiles  were  obtained  by  hot-wire  traverse  5  diameters  upstream 
of  the  bend  and  at  90°. 

The  remaining  tests  were  made  with  70  diameters  of  inlet 
flow  development  in  order  to  establish  essentially  fully-developed 
conditions  at  entry  to  the  bend.  Limited  temperature  and 
velocity-field  data  were  taken  at  the  same  Reynolds  number  as  the 
earlier  test;  these  indicated  no  significant  variations  in  flow 
structure.  In  the  limited  time  available  attention  was  mainly 
directed  at  obtaining  data  on  heat  transfer  coefficients  over  as 
wide  a  range  of  Reynolds  number  as  possible  rather  than  focusing 
on  a  very  detailed  mapping  of  the  interior  temperatures  at  one 
Reynolds  number.  Data  were  thus  obtained  at  nominal  Reynolds 
numbers  of  9  x  10 3  and  9  x  104. 


2.1.3  The  Experimental  Results  Figure  1  compares  Berkeley  and 
UMIST  measurements  of  streamwise  velocity  at  90°  around  the  bend. 
A  striking  feature  of  the  Berkeley  measurements  was  the  double 
peak  in  streamwise  velocity  (a  feature  that  neither  their  nor  our 
computations  reproduce).  The  velocity  data  obtained  at  UMIST 
confirm  this  feature  of  the  Berkeley  experiments.  The  small 
differences  between  the  two  sets  of  readings  could  well  be  due  to 
a  positioning  error  in  the  UMIST  data.  (Our  original  intention 
had  been  to  take  only  a  single  bottom-to-top  traverse  along  the 
mean  radius  of  curvature  but,  to  allow  more  direct  comparison  with 
the  Berkeley  data,  which  were  obtained  from  slde-to-side 
traverses,  rather  rudimentary  adaptations  were  made  to  allow 
traverses  at  other  radii).  It  was  concluded  therefore  that  we 
could  assume  that  the  flow-field  in  the  UMIST  apparatus  was  indeed 
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essentially  the  same  as  documented  in  the  geometrically  similar 
half-scale  Berkeley  apparatus. 

O  A 

Temperature  profiles  across  the  duct  are  shown  for  45  ,  90 
135°  and  190°  in  figure  2.  These  mimic  closely  the  corresponding 
velocity  profiles.  For  example,  by  45°  the  "peak"  temperature 
(actually  the  coolest  point  in  the  stream)  has  shifted  right  of 
centre,  while  at  90°  there  is  a  deep  trough  in  the  temperature  in 
just  the  same  position  as  the  trough  in  velocity  shown  in  figure 
1.  The  remains  of  these  troughs  are  still  clearly  present  at  l'S° 
and  are  faintly  visible  even  at  180°. 

The  development  of  the  Nusselt  number  around  the  bend  is 
shown  in  figure  3.  Up  to  entry  to  the  180°  section  the  Nusselt 
number  is  nearly  uniform  along  all  walls.  By  45°,  however,  the 
level  of  Nu  on  the  inner  line  of  symmetry  has  fallen  by  25%  while 
the  average  level  has  risen  by  about  10%.  By  90°  there  is  more 
than  a  2:1  ratio  in  the  Nusselt  numbers  recorded  on  the  inner  and 
outer  symmetry  planes.  This  ratio  remains  nearly  constant  over 
the  second  half  of  the  bend  though  the  absolute  level  of  Nu 
diminishes  slowly.  There  remains  a  significant  difference  in  heat 
transfer  coefficient  even  10  diameters  downstream  of  the  bend,  the 
heat  transfer  coefficient  on  the  outer  bend  wall  being  30-50% 
higher  than  on  the  inner  wall.  Throughout,  the  heat  transfer 
coefficient  on  the  side  wall  follows  closely  that  on  the  outer 
wall . 


2. 2  Computational  Program 


2.2.1  Numerical  and  physical  Model  The  computer  program 
embodying  a  three-dimensional  semi-elliptic  solving  scheme  for  the 
averaged  equations  of  motion  was  provided  by  Professor  Humphrey. 
It  is  based  on  the  discretizational  and  programming  strategy  of 
the  TEACH  family  of  computer  codes  save  that  convective  transport 
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is  h*re  represented  by  quadratic  upstream  weighting  [  "f  ]  .  The 
code  incorporates  as  alternatives  the  k~e  Boussinesq  viscosity 
and  the  algebraic  stress  models  of  turbulence,  both  of  which 
utilize  the  standard  transport  equations  for  the  turbulence  energy 
and  its  dissipation  rate  e. 

The  effort  at  tWIST  initially  consisted  of  a  very  thorough 
re-checking  of  the  algebraic  and  coded  forms  of  the  transport 
equations  for  the  flow-field  variables.  This  exercise  proved  very 
useful  in  establishing  confidence  that  what  was  inevitably  a  large 
and  complex  code  was  free  from  coding  error. 

Thereafter  the  computer  program  was  adapted  to  incorporate 
a  new  wall-function  treatment  for  the  viscosity-affected  zone 
between  the  near-wall  node  and  the  wall  itself  [ft],  A  further 
adaptation  was  the  inclusion  of  the  streamline-curvature 
modification  to  the  transport  equation  fore[9]  proposed  on  the 
basis  of  two-dimensional  flows.  Further  details  are  given  in  ref. 
[1] appended.  A  solving  routine  was  also  included  for  the  thermal 
energy  equation.  Since  the  fluid  properties  were  taken  as 
independent  of  temperature,  this  equation  was  solved  only  after  a 
converged  solution  had  been  obtained  to  the  velocity  field. 

7.2.2  Comparison  of  Computed  and  Measured  Behaviour  Initial 
computations  were  made  using  a  70  x  J7  non-uniform  grid  to  map  the 
cross-section  of  the  duct  with  110  streamwise  planes,  50  of  which 
were  located  on  the  bend  itself  (i.e.  a  1 0  spacing).  The  nodal 
density  was  subsequently  refined  to  15  x  7°  in  the  cross-section 
while  over  the  first  45°  of  the  bend  computational  planes  were 
spaced  at  1$°  Intervals.  These  refinements,  while  leading  to 
changes  up  to  a  maximum  of  10%  in  the  secondary  flow,  had  only  a 
very  weak  effect  on  the  streamwise  velocity  profiles. 

The  situation  regarding  computations  can  most  effectively 
be  conveyed  by  reference  to  the  90°  plane  (ref.[l]  in  the  Appendix 
shows  further  details).  Prom  figure  i  it  is  evident  that  the 
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large  trough  in  the  streamwise  velocity  shown  by  the  experiments 
is  not  reproduced  by  the  computations.  There  are  unfortunately  no 
data  between  the  43°  and  90°  stations  but  the  computations  seem  to 
indicate  that  the  trough  has  arisen  through  the  secondary  velocity 
separating  off  the  inside  wall  somewhere  between  these  two 
stations  thus  giving  a  secondary  recirculation,  figure  5,  in  the 
opposite  sense  than  usual.  It  must  be  said  however  that 
measurements  at  90°  show  no  sign  of  such  a  reversed  eddy  so  this 
explanation  remains  tentative.  The  behaviour  shown  in  figure  4  is 
nearly  the  same  as  that  obtained  at  Berkeley  by  Chang  [19]  using, 
in  both  cases,  the  k-*e  Boussinesq  viscosity  model.  The  relative 
insensitivity  of  the  results  to  grid  refinement  and  to  the  choice 
of  wall  function  suggests  that  the  poor  predictions  are  mainly  due 
to  the  turbulence  model  used  in  the  main  flow  region.  Yet  it  is 
found  [1]  that  very  little  effect  resulted  from  introducing  the 
curvature  correction  to  the  transport  equation  or  from  adopting 
the  algebraic  stress  model.  Although  Mr.  Johnson  has  now 
concluded  his  research.  Professor  Young  Don  Choi,  an  academic 
visitor  mentioned  in  the  footnote,  continues  .to  address  the 
problem. 


3  ROUND  SECTIONED  DUCT 


3.1  Experimental  Program 

The  experimental  work  on  the  round  sectioned  duct  has 
greatly  benefited  from  extensive  interaction  with  Professor  J.W. 
Baughn  of  the  Mechanical  Engineering  Department,  University  of 
California  Davis.  We  concluded  that  the  attainment  of  a 
uniform-heat-flux  boundary  condition  was  impracticable  for  this 
particular  geometry  and  so  the  apparatus  design  has  been  based  on 
a  uniform  wall  temperature  rig.  The  ]*n°  bend  section  has  been 
fabricated  in  two  halves  (figure  machined  from  solid  blocks  of 
aluminium.  The  ratio  of  bend  radiustpipe  diameter  is  the  same  as 


*  Current  work  at  UKIST  by  Professor  Young  Don  Choi  with  twice  as  many  cross- 
sectional  nodes  and  50%  more  streamvise  planes  has  led  to  only  minor  improve¬ 
ments. 
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for  the  square  sectioned  duct,  3.35:1.  Inlet  and  outlet  tangent 
sections  each  *0  diameters  in  length  made  of  high-grade  aluminium 
tubing  are  installed.  At  the  downstream  end  of  this  assembly, 
connection  is  made  to  the  flow-metering  section  and  fan  exhaust 
used  for  the  square  duct  experiment. 

Heat  transfer  rates  to  the  duct  wall  are  to  be  obtained  by 
heat  flux  meters  based  on  the  design  of  Professor  Baughn  and  one 
of  his  students  at  UC  Davis,  ref.  [11]  .  The  holes  into  which  the 
meters  will  be  cemented  are  clearly  visible  in  figure  5.  A  sketch 
of  the  meter  is  shown  in  figure  7.  The  principle  of  operation  is 
that  heat  is  supplied  via  resistors  to  the  flat  metal  cone  at  a 
rate  that  is  just  sufficient  to  keep  the  cone  temperature  the  same 
as  that  of  the  pipe.  At  present  a  prototype  meter  has  been  built 
and  is  undergoing  testing. 


3.2  Computational  Program 

3.7.1  Numerical  Solution  Procedure  In  the  19PI/R7  Annual  Report 
we  reported  that  a  semi -el 1 i pt ic  procedure  for  flow  through 
toroidal  ducts  had  been  developed  and  showed  its  successful 
application  to  one  of  the  laminar-  flow  cases  of  Agrawal,  Talbot 
and  Gong  [12]  at  a  Dean  number  of  1P3.  During  19P2/93  the 
procedure  has  been  applied  to  duct  flows  at  successively  higher 
Dean  numbers  and  this  has  pointed  the  need  for  further 
refinements. 

An  Important  step  has  been  the  replacement  of  the  SIMPLE 
algorithm  by  the  more  recent  SIMPLER  procedu-re,  Patankar  [17]  The 
former  approach,  while  successful  enough  when  curvature  terms  were 
moderate,  did  not  succeed  in  procuring  convergence  as  the  Dean 
number  was  successively  raised.  Details  of  how  SIMPLER  has  been 
implemented  in  the  present  semi-elliptic  scheme  a  re  given  in 
reference  [7]  and  [3]  contained  ijj  the  Appendix. 
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A  numerical  refinement  of  a  different  kind  has  been 
prompted  by  the  problem  of  computing  turbulent  flow*  Because  of 
the  very  strong  secondary  flow  that  is  generated  close  to  the 
wall,  it  was  felt  strongly  desirable  to  discard  the  habitually 
used  wall-function  approach  and  instead  carry  out  the  integration 
to  the  wall  itself.  However,  there  was  insufficient  core 
available  to  us  to  add  the  further  ten  or  so  radial  lines  needed 
to  resolve  the  viscous  and  buffer  regions +  .  The  problem  was 
removed,  however,  by  recognizing  that  across  the  region  where 
viscous  effects  were  significant,  the  pressure  variation  would  be 
adequately  obtained  from  assuming  radial  equilibrium.  In  a  semi- 
elliptic  solver  it  is  the  pressure  that  limits  the  mesh  density 
because  only  this  variable  requires  three-dimensional  storage. 
Thus,  we  put  in  additional  velocity  nodes  across  the  near-wall 
sublayer  but  no  corresponding  nodes  for  the  pressure.  Further 
details  may  be  found  in  [4]  in  the  Appendix. 

The  computations  for  turbulent  flow  have  so  far  adopted  the 
standard  k—e  Boussinesq  viscosity  model  in  the  fully  turbulent 
region  matched  to  the  Van  Driest  form  of  the  mixing  length 
hypothesis  across  approximately  the  51  of  the  flow  nearest  to  the 
wall.  Further  details  are  given  in  ref.  [3]  in  the  Appendix.  In 
solving  the  thermal  energy  equation  to  predict  heat  transfer 
coefficients,  a  uniform  turbulent  Prandtl  number  of  n.o  is  adopted 
throughout  the  pipe. 

3.2.2  Examples  of  Applications  Although  the  main  interest  in 
the  project  is  in  turbulent  flow,  a  thorough  testing  of  the 
computational  scheme  for  laminar  flow  was  felt  desirable  in  view 
of  the  several  sets  of  data  available.  In  continuation  (from  our 
19H1/A2  report)  of  the  study  of  the  Agrawal  data  [12],  figure  R 
compares  streamwise  velocity  profiles  at  two  sections  for  a  Dean 
number  of  553.  The  computations  started  at  the  entry  to  the  bend 
with  a  uniform  inlet  velocity  assumed.  The  grid  employed  was 
20  x  20  x  100.  Close  agreement  is  displayed  between  the  computed 
and  measured  behaviour.  A  further  case  examined  was  that  of 


In  the  present  work  this  region  is  spanned  with  the  mixing  length  hypothesis. 
More  then  10  nodes  would  probably  be  required  if  a  more  elaborate  closure 
were  adopted. 
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Enayet  et  al  [14]  in  a  90°  bend  with  a  radius :diameter  ratio  of 
only  2.9:1  giving  a  Dean  number  of  464.  in  this  case  the  bend  was 
preceded  by  a  straight  entry  section  and,  as  a  result,  at  entry 
the  boundary  layer  thickness  was  about  0.5  times  the  pipe  radius. 
The  higher  Dean  number  and  the  presence  of  relatively  thick 
boundary  layers  leads  to  a  stronger  secondary  flow  being 
established.  Figure  9  indicates,  however,  that  a  satisfactory 
numerical  simulation  is  nevertheless  obtained. 

Enayet  et  al  [14]  also  measured  the  development  of 
turbulent  flow  through  the  same  duct  and  our  computations  of 
streamwise  velocities  in  this  flow  are  shown  in  figure  10. 
Agreement  is  now  less  complete  than  for  the  laminar  flow  case  and, 
in  particular,  the  trough  in  the  velocity  profile  at  5n°  and  75° 
is  only  qualitatively  predicted.  Nevertheless,  the  measured 
behaviour  is  simulated  far  more  satisfactorily  than  for  the  square 
sectioned  duct  discussed  in  Section  2. 

Computations  of  the  heat  transfer  behaviour  are  provided  in 
ref.  4  in  the  Appendix.  Comparisons  with  the  experiments  of 
Seban  and  McLaughlin  [15]  of  fully-developed  flow  in  a  coil 
suggest  that  our  scheme  predicts  accurately  the  rise  in  the 
circumferentially  averaged  heat  transfer  coefficient  but  that  the 
augmentation  on  the  outside  of  the  bend  is  underestimated  as  is 
likewise  the  damping  around  the  inside.  This  is  qualitatively  the 
type  of  behaviour  to  be  expected  with  a  Boussinesq-type  model 
since  it  does  not  mimic  the  great  sensitivity  to  streamline 
curvature  that  real  turbulence  displays.  The  thermal  field 
computations  of  the  2.9:1  90°  bend  of  Enayet  et  al  [14]  indicate  a 
five-fold  variation  of  local  heat  transfer  coefficient  around  the 
bend  at  75°.  The  experiment  was  purely  concerned  with  the  flow 
field  -  no  heat  transfer  was  involved  -  but,  on  the  basis  of  the 
comparison  with  the  Seban-McLaughlin  data,  it  would  seem  probable 
that  in  that  geometry  as  mu  as  a  ten-fold  variation  in  heat 
transfer  coefficient  may  actually  occur  between  the  inside  and 
outside  of  the.  bend. 


Future  Work 
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4.1  Experimental  Research 

The  principal  experimental  work  concerns  the  commissioning 
of  the  ci rcular-sectioned  tube  apparatus.  First  tests  on  the 
assembled  rig  will  begin  during  April  19P4.  The  initial  testing 
will  be  concerned  with  the  in  situ  performance  of  the  heat  flux 
meters,  tested  in  the  straight  entry  section  against  well- 
established  data  of  developing  flow  in  a  circular  tube.  The 
pattern  of  experiments  will  be  similar  to  that  in  the  square- 
sectioned  duct:  a  detailed  set  of  measurements  of  wall  flux  and  a 
mapping  of  the  temperature  field  at  the  same  Reynolds  number  as 
the  Berkeley  circular  tube  experiments.  Some  limited  confirmatory 
velocity  field  data  will  also  be  gathered.  Thereafter  wall  flux 
data  will  be  obtained  at  the  extremes  of  Reynolds  number 
accessible  to  the  fan. 

Looking  further  ahead  (further  than  the  period  covered  by 
the  present  grant},  we  will  aim  to  obtain  from  the  two  l«n°  bend 
apparatuses  several  additional  sets  of  experimental  data.  These 
will  include  explorations  where  the  flow  entering  the  bend  is 
essentially  at  uniform  velocity  and  studies  at  the  low-Reynolds- 
n umber  end  of  the  turbulent  regime  where  partial  laminarization 
may  occur  in  passage  around  the  bend. 


4.2  Computational  Work 

Our  present  efforts  are  directed  at  incorporating  an 
algebraic  stress  model  into  the  circular  bend  code.  There  is 
reasonable  optimism  that  useful  improvements  in  the  fidelity  of 
the  computer  simulation  will  result.  The  lingering  doubt  concerns 
the  fact  that,  when  included  in  the  square  duct  code,  the 
algebraic  stress  model  did  not  produce  better  results.  Against 
this  discouraging  fact,  however,  may  be  set  many  striking 
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successes  (recorded  in  the  literature  of  the  last  ten  years)  from 
applying  algebraic  stress  models  to  curved  flows.  Moreover,  in 
contrast  to  the  fairly  satisfactory  prediction  of  the 

O 

round-sectioned  90  bend  with  the  k~e  Boussinesq  viscosity  model, 
the  square-sectioned  duct  predictions  with  the  same  turbulence 
model  go  completely  wrong  beyond  49°  of  arc.  Our  current  view  is 
that  the  poor  accuracy  of  predictions  in  the  square  sectioned  duct 
is  mainly  due  to  weaknesses  in  the  dissipation  rate  equation. 
Over  the  remaining  period  of  the  grant  Professor  Choi's  efforts 
will  be  directed  at  improving  this  aspect  of  the  modelling. 
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Figure  1  Comparison  of  normalised  mean  streanwise  velocity  et 
90°  in  the  bend;  O  "  date  Cheng  et  el  (1983) ,  31  Du 
inlet.  Re  ■  56,690;  +  -  date  present  study,  30  0- 
inlet,  Re  ■  57,300;  •  -  data  present  atudy,  72  D_ 
inlet.  Re  ■  56,500. 
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Figure  2  Turbulent  flow  measurement*  of  seen  temperature  at 

0  »  45  ,  135  and  180°  in  the  band  ,  non-dimensional ised 
with  reapaet  to  tha  wall  temperature  (T„)  on  the  outer 
wall  at  2x/D_  ■  1  and  the  bulk  temperature  (T-)  for 
Re  •  55,800,  inlet  tangent  30 Dg 
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Figure  3(a)  Nusselt  number  computed  from  experiment  at  5 

upstream  and  0*  and  45*  in  the  bend;  o  -  outer 
wall,  a  -  bottom  wall,  *  -  inner  wall;  inlet 
tangent  30  DH;  Re  *  56,030. 
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TURBULENT  MOMENTUM  AND  SEAT  TRANSPORT  IN 
SUM  THROUGH  A  ISO0  BEND  OF  SQUARE  CROSS-SECTION 
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ABSTRACT 

Tha  paper  reports  flow  aad  haat  trangfer  predictions 
of  turbulent  flow  in  paasaga  around  a  180°  square- 
sectioned  band.  Tha  numerical  raaulta  ara  obcainad  from 
a  finite-volume  discretization  of  cha  aani-allipcic  fora 
of  cha  Naviar  Stokaa  and  energy  aquaciona.  Tha  turbulant 
acraaaaa  ara  rapraaancad  by  cha  k-t  Bouaainaaq  viscosity 
aodal  boch  in  ica  acandard  fora  and  wich  a  acraaalina 
curvature  cor race ion.  Sarioua  diffarancaa  baewaan 
experiment  and  pradiccion  axiat  for  boch  forms  and 
auggaaeiona  ara  aada  for  chair  origin.  Haac  cranafar 
computations  indicaca  chac  froa  43°  -  180°  cha  aacondary 
flow  provokes  ac  laasc  a  2:1  eircuafarantial  variacion 
in  haac  cranafar  coafficianc  around  cha  duct  pariaacar 
and  chac  diffarancaa  of  40Z  baewaan  cha  onan  haac 
cranafar  coafficianc  in  aach  wall  paraiac  at  laaac  10 
diaamters  downstream. 


1  nSTRODUCTIO:: 

Ca^ucara  and  flew  aolucion  schemes  hava  davalopad 
co  a  poine  whara  aarioua  nuaarical  acudiaa  of  coovaceiva 
haat  cranafar  in  complex,  3-dimansional  flowa  can  now  ba 
■ada.  But,  for  turbulant  flow,  what  1ml  of  accuracy 
can  ona  expect  froa  auch  a  simulation?  Tha  quaacioc  ia 
of  groat  practical  interest  for  if  cha  accuracy  can  ba 
raliad  on  wich in  cha  allowable  colarancaa  cha  coat  of 
computation  -  typically  a  f aw  hundred  dollars  par  run 
-  will  oftan  bo  nagligibla  coaparad  with  the  coat  of 
getting  tha  inforaatioo  froa  experiment.  This  ia  cha 
question  our  rasaarch  on  flow  around  180°  bands  has  baan 
addressing. 

Tha  !80°-bend  flow  has  several  qualities  that  aaka 
it  wall  suited  as  a  bench-mark  cast  ease.  It  has  vary 
strong  practical  connections,  especially  in  haat 
exchangers,  yet  its  topography  ia  relatively  simple. 

This  simplicity  naans  chac  cha  flow  boundary  conditions 
can  ba  easily  and  unambiguously  raproducad  by  a  comeucor; 
it  also  naans  that  obtaining  tolerably  accurate  numrioal. 
solutions  (whatever  nay  ba  the  frailties  of  tha  physics) 
is  a  target  within  sights  -  though  Chare  nay  ba  argtasnt 
about  whacker  ie  is  yac  within  range.  It  is  a  wore 
challenging  flow  than  tha  90°  band  chat  provided  one  of 
Cha  case  cases  at  the  1981  Stanford  Confaroncs  (Elina 
at  al,  1981,  1982)  because  turbulent  stresses  gonarated 
by  tha  strong  secondary  flow  hava  longer  to  act  on  cha 
■tan  flow.  Moreover ,  ad  nest  import  ant  ly,  detailed 
esporimancal  data  era  available  (Chang  at  al,  1983)  with 
which  to  draw  comparison. 

Computations  of  flow  around  tha  90°  tquara-saccioned 
band  adopted  for  tha  Stanford  Conference  hava  boon 
reported  by  Humphrey  at  al  (1981,  McDonald  (1982), 
Abddlaaguid  ac  al  (19(2),  Radi  at  al  (1982),  Moora  tad 
Moots  (X9S2)  and  Chang  ac  si  (1982);  tha  first  three 
employ  a  discretisation  of  tha  full  Naviar  Stokes 
aquations,  cha  last  two  adept  a  aeni -elliptic  formulation 
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in  which  only  tha  pressure  field  is  stored  over  cha  full 
domain  (Pratap  and  Spalding,  1975)  and,  necessarily, 
sereamriae  derivatives  in  tha  momentum  equations  are 
discarded.  All  but  McDonald's  and  Moore's  studies 
•■ployed  cha  same  mathematical  Mel  for  the  turbulent 
stress  field  -  tha  k-c  Boussineaq  viscosity  modal.  Tat 
the  flow  patterns  predicted  by  these  schemas  bear  only  a 
qualitative  rasaablanca  to  each  ochar  or,  for  that 
matter,  co  tha  experimental  data.  The  fully  three- 
dimensional  discretizations  hava  of  necessity  co  use  a 
coarse  grid  which  has  severely  limited  tha  nuamricel 
accuracy  available.  The  semi-ellipcic  schemas  which  have 
allowed  80  or  more  streamwise  planes  (not  all,  of 
course,  in  che  bend  itself)  achieve  somewhat  better 
overall  agreement  with  data. 

The  computations  of  Moore  and  Moore  (1983)  and 
McDonald  (1983)  extend  through  the  viscous  sublayer  to 
the  well  allowing  a  better  numerical  resolution  of  the 
near-wall  region  than  tha  k-t  treatments  (which  apply 
"wall  functions"  to  bridge  che  5Z  of  the  flow  nearest 
the  wall).  This  fine- grid  approach  was  also  followed  in 
tha  work  of  Couataix  at  al  (1983)  for  tha  same  conference 
though  they  make  only  a  two-dimensional  inviecid 
calculation  for  tha  pressure  field.  Tha  maeerical 
simplification  this  brings  is  considers^’  -i—'t  a  three- 
dimensional  marching  scheme  may  ba  adopt  d;  iu  view  of 
tha  poor  agreement  obtained,  however,  ft  remains 
questionable  whether  this  basis  for  obtaining  the  pressure 
is  useful  in  bands  of  small  aspect  ratio  with  substantial 
curvature.  Of  course,  more  nodes  near  the  wall  means 
fewer  elsewhere;  moreover,  instead  of  calculating  the 
dissipation  rata  t  from  a  transport  aquation  these 
groups  obtained  ic  via  a  prescribed  length  scale  distrib¬ 
ution  -  with  much  uncertainty  as  to  the  appropriate 
proscription. 

In  tbe  intervening  two  years  since  computations  for 
tha  Stanford  Conference  ware  made  the  authors  have  cont¬ 
inued  to  giva  attention  to  curved  ducts  but,  for  reasons 
given  earliar,  to  che  180°  bend  case. 

A  detailed  mapping  of  the  velocity  field  by  laser 
anenemetry  has  been  made  at  Berkeley  for  the  curved  duct 
shown  in  figure  1.  These  data  aad  preliminary  ntmerical 
computations  ara  to  appear  in  Chang  at  si  (1983).  Tha 
present  contribution  provides  s  more  refined  set  of 
eompueacions  of  this  flow  and  examines  the  influence  of 
wall  boundary  conditions  and  turbulence  modal  on  the 
computed  flow  pattern.  It  alto  reports  solutions  of  the 
enthalpy  equation  thus  providing  predictions  of  the 
circumferential  distribution  of  Nussolt  lumber  around 
the  duct  perimeter. 

2.  SUMMIT  OF  NUWR1CAL  SCREW 

Describing  Differential  (auctions 

The  stationary,  turbulent,  ineo^ressible  flow  of 
fluid  through  a  curved  duct  of  constant  rectangular 
cross-seetion  is  conveniently  described  through  conserv¬ 
ation  equations  in  cylindrical  coordinates.  With 
coordinates  x  sad  r  mapping  the  duct  cross-section  and 


9  the  sngls  of  progress  along  cha  dues,  tba  describing 
Man  flow  aquation*  aav  ba  vricten: 


00.  ♦  VC.  *  (W/r)U.  --ip.*  (2u0  -  ?). 
x  r  vox  x  x 

♦  r"1(r(v(0r  ♦  Vx)  -  uv)f  (l) 

r-nomsntun 

0VX  ♦  Wr  -  H*/r  ♦  (W/r)V9  -  -  i  Pr  ♦  (v(0r  ♦  Vx)-OT)x 
♦  r_1(r(2vVr  -  7f)r  -  2vV/ r*  (2) 


9  nomsntun 

OH  *  VH  *  VH/r  ♦  (H/t)H  -  -  i  P,  *  (v(H  ♦r"10.)^«). 
x  r  o  or  o  xox 

♦  r'2(ri(v(r(H/t)r  *  r_1V9)  -  w)>r  (3) 

continuity 

V  *  V/r  ♦  r'V  *0-0  (4) 

r  9  * 

enthalpy 

UT  ♦  VT  *  (U/ r)T  •  ((v/a)T  -  u?) 
x  r  o  xx 

♦  r*l(r((u/q)rr- 5?))r  (3) 


Hara  U,V  and  U  arc  qta  naan  velocity  components  in  the 
x,  r,  9  directions,  u;  uv,  ece.  are  the  correspondingly 
defined  Reynolds  scresaea,  T  is  the  naan  eaaperstura  and 
t  the  temperature  fluctuation.  The  quantities  v,  s  and 
o  are  respectively  the  fluid  kinematic  viscosity, 
density  and  Prandtl  nuafaer.  The  subscripts  x,  r  and  9 
denote  partial  differentiation  with  respect  to  tba  space 
coordinate  in  question. 

The  turbulent  stresses  are  obtained  fron  the 
Boussinesq  stress-strain  formula  which,  in  cylindrical 
coordinates,  i^ilies: 

-  u2  •  2vt  Dx  ”  •jk;  "  2vy.  vr  “  ■jk 

-  «*  •  vt(Hx  ♦  r-10e) }  -v*  -  vT(r(W/r)r  ♦  r”1"^) 

-  uv  •  vT(0r  ♦  Vx)  (6) 

where,  in  accordance  with  cha  k-c  nodal 

VT  “  cuki/e  •  <7) 

k  being  the  curbulenca  energy  and  c  its  dissipation 
race.  These  quantities  are  chenaelves  found  fron  trans¬ 
port  aquations  solved  simultaneously  with  the  naan  flow 
variables: 

0kx  *  Vkr  *  (H/r)k#  -  r_l(r(  V>r  *  (  kx>x 

♦  v,(2(0*  *  v£  ♦  Tm2v*  *  o  V  *  r_10.H_  *  r’Vvj 

I  x  y  o  r  X  9  X  vr 

♦  o*  ♦  r”*0?  *  v*  ♦  r m*v'i  *  v'i  ♦  h* 

f  V  |  V  I  t 

♦  *  V)  -  t‘hnnt  *  2t«r  -  U))  -  t  (8) 


0ax  ♦  Rtf  *  Wr>*t  "  V^T 
*  eclvT  (2(0*  *  0*  *  r-V#  *  0^  * 


*r»r*(  V, 

*'ty,  * 


♦  c|  *  r_20|  *  V*  ♦  r“*7 j  *  W*  ♦  V* 

*  2r_2V  (2He  ♦  V)  -  r~h/(2Ve  *  2rHf-  H))  -  et2e*/k 

(9) 

The  a^irical  coefficients  are  assigned  the  usual  values 
optimised  by  Launder  et  al  (1973) : 

eu  •  0.09;  c  •  1.44;  C(J  •  1.92;  ^  •  1.0;  o(  •  1.22 

The  aemi-alliptic  truncation  has  been  applied  to  all 
Chase  aquations  in  that  terns  containing  second  deriv¬ 
atives  on  9  are  dropped. 

The  Oifferenca  Equations  and  their  Solutions 

Fiaice-dliferance  forms  of  the  transport  equations 
wars  derived  by  integration  oyer  discretized  volumes  in 
the  flow  domain  following  broadly  Pracap 1 s  (1975) 
guida lines.  As  is  now  cus tonary,  the  velocity  conponents 
and  pressure  are  stored  on  a  staggered  mesh .  Streanwise 
convective  transport  is  approxiaacad  by  upwind  differences. 
In  cha  cross-sectional  plana,  however,  both  upwind  and 
quadratic  upwind  (QUICK)  options  are  included.  The 
latter  seheas  devised  by  Leonard  (1979)  has  been  tested 
by  Han  at  al  (1981)  for  turbulent  flow  in  a  driven 
cavity  and  found  to  be  distinctly  better  than  an  upwind 
approximation.  (In  fact,  so  far  as  the  cross-streaa 
conponents  are  concerned,  the  flow  around  a  bend  is  very 
like  a  driven  cavity).  The  nethod  of  iaplenenting  Che 
schema  in  the  mane ri cal  algorithm  is  as  given  by  Han  et 
al  (1981). 

In  ehe  course  of  iteration,  the  adjustments  to  the 
pressure  field  in  response  to  suss  imbalances  for  the 
control  volumes  surrounding  eech  pressure  node  is 
essentially  as  proposed  in  the  SIMPLE  algorithm  (Patankar) 
and  Spalding,  1972)  axeape  chat,  with  the  velocity  field 
held  on  only  two  successive  9  planes,  reorganization  is 
required  (Pracep ,  1975).  At  any  x-y  plane  the  U  and  V 
velocity  conponents  are  solved  first.  The  streawise 
momentum  equation  is  solved  next  to  obtain  the  H  velocity 
(displaced  half  a  cell  down-stream  aa  a  result  of  the 
'staggering')  using  new  values  of  v  and  U  in  the  convect¬ 
ive  tense;  finally,  perturbations  to  the  pressure  field 
are  introduced  in  conjunction  with  re-sdjustaants  to 
the  currant-plane  U-V  field.  This  procedure  is  applied 
at  all  planes  beginning  at  the  upstream  boundary  and 
stopping  downstream,  succassively  overwriting  'upstream' 
velocities  by  currant  values.  On  completing  such  a 
pass  over  the  domain  the  confutation  starts  over  again 
at  the  upstream  and  unless  the  residual  error  is  snail 
enough  that  convergence  is  signalled. 

Because  the  confutation  as  outlined  above  has  had 
to  naka  extensive  use  of  upstream  values  (rather  than 
currant  plans  values)  of  velocity  in  evaluating  convection 
coefficients  and  sources,  a  certain  upscream  bias  is 
introduced  into  the  solution  if  a  purely  narching  treat- 
neat  is  followed.  For  this  reason,  as  the  cogitation 
approaches  its  apparent  solution,  it  is  necessary  to 
introduce  iteration  on  che  velocity  components  at  each 
step.  That  is  to  say,  whan  current  plane  values  have 
been  obtained  the  aquations  ere  resolved  using  current- 
values  aa  appropriata  in  -forming  coefficients  and 
source  terms.  Approximately  45  passes  were  needed, 
starting  from  an  assumed  uniform  pressure  field  to  obtain 
converged  results;  this  wss  deemed  to  have  been  achieved 
whan  cha  magnitude  of  the  mass  errors  siaaeed  over  every 
cell  of  the  domain  fell  below  IX  of  the  entering  aass 
flow.  (Other  studies  have  typically  accepted  ness 
errors  of  0.5X  par  plena  which  is  larger  by  e  factor  of 
SO  than  tolerated  here). 

The  flow  field  generated  in  the  duet  is  symmetric 
about  the  mid-plane  of  the  cross-section  lying  in  the 
plane  of  the  bend.  Computations  were  thus  extended 
over  just  one  half  of  the  duct,  ehe  cross-section  being 
napped  by  a  15x23  interior  grid  for  nest  of  the  results 
presented  here-under  with  the  nesh  expanding  nildly  from 
each  well.  The  conputstions  began  seven  hydraulic 
diMstars  upstream  of  eha  bend  sad  extended  11. S 
downstreams  ehis  region  was  covered  by  a  eotel  of  117 
streaMrlaa  planes. 
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boundary  Condition* 

Along  the  sySstry  plana  th*  (radiant i  of  all  but 
ona  of  tha  dapandant  variablaa  wara  aat  to  aaro;  tha 
valua  of  U,  tha  velocity  normal  to  this  plana,  uai  sada 
aaro.  On  tha  thraa  aides  boundad  by  tha  duct  vail,  vail 
functiona  vara  employed  to  supply  appropriata  near-well 
aourcaa  and  sinks  to  tha  various  dapandant  variablaa, 
Basidas  tha  'standard'  vail  traacaant  habitually  alloyed 
in  codas  davalopad  by  thoaa  aaaociatad  vith  tha  lm>*rial 
Collaga  school  (Launder  and  Spalding,  1974)  a  more 
alaborata  varsion  originating  in  tha  prasant  vork  haj 
baan  usad.  It  is  an  extension  of  tha  schaaas  of  Chieng 
and  Laundar  (1980)  and  Johnson  and  Laundar  (1982).  In 
ralacion  to  tha  prasant  study  its  aost  significant 
feature  is  that  tha  vail  friction  opposing  tha  secondary 
taotion  is  obtained  independently  of  tha  streaarise  val- 
ocity  by  nerforming  tha  integration  V  -  /(t^ )dx 

batvaen  tha  vail  and  tha  first  node.  By  contrast,  tha 
IC  vail  treatment  assumes  that  tha  resultant  naar-vall 
velocity  parallel  to  tha  wall  obeys  tha  usual  logaritbtdc 
lav  (tha  streaatrise  and  cross-stress  vail  stress  cov¬ 
enants  are  then  obtained  by  resolving  appropriately). 

Flow  inlet  conditions  are  detailed  in  tha  next 
section. 


3.  COHPPIATiniS  BMP  BTOMB  CagABEP 

Tha  experimental  data  providing  the  basis  for  this 
cosparison  are  thoaa  reported  by  Chang  at  al  (1983)  for 
a  duct  in  which  tha  naan  radius  of  tha  band  was  3.33 
hydraulic  disasters.  Tha  flow  entering  tha  band  had 
davalopad  through  a  straight  entry  section  of  31  D„  after 
being  passed  through  a  aariea  of  screens  to  promote 
shear  layer  development.  Thus,  while  tha  flow  had  not 
become  fully  davalopad  there  was  no  inviscid  core 
remaining  whan  tha  flow  encountered  tha  band.  Tha 
experiments  were  taken  at  a  bulk  Reynolds  number  of  56700. 

A  parallel  experiment  at  UMIST  is  underway  which 
reproduces ,  so  far  as  we  are  able  to,  tha  Berkeley  teat 
conditions.  Tha  apparatus  disenaions  are  twica  thoaa  of 
tha  Berkeley  rig  and  since  air  rather  than  water  is 
employed,  velocities  need  to  be  increased  by  a  factor  Of 
8  to  maintain  the  same  laynolda  number.  While  the  main 
output  from  tha  OMIST  study  will  ba  convective  beat 
transfer  data  it  has  served  to  provide  checks  on  an 
unexpected  feature  of  the  Berkeley  maaaurawnts.  At  90° 
around  the  bend  thair  atraamviaa  velocity  profile  along 
a  radial  line  had  exhibited  a  pronounced  double 
maximum.  This  feature  while  vat  strongly  prasant  on 
the  plant  of  symmetry  was  still  evident  along  tha  line 
midway  between  tha  symmetry  plana  and  the  and  wall. 

Figure  2  coheres  the  laser  anemometer  profiles  reported 
by  Chang  at  al  (1983)  with  those  obtained  vith  a  pair  of 
slant  hot  wires  at  UKIST.  It  is  seen  that  tha  double¬ 
peak  faatura  it  present  in  both  eats  of  data;  indeed, 
there  is  vary  satisfactory  agreement  batvaen  tha  two 
realisations  of  this  flow. 

Computations  were  started,  as  noted,  seven 
hydraulic  diameters  upstream  of  the  bend  using,  as 
initial  conditions  tha  velocity  and  turbulence  energy 
profiles  of  Mailing  and  White lm.  (1976).  The. turbulence 
energy  dissipation  rata  was  assigned  as:  t  »k3,2/i  where 
the.length  seals  was  assigned  as  tbs  smaller  of 
*/e  /4  times  tha  distance  to  tha  nearest  vail  or  .375D_ 

Tha  former  it  consistent  with  a  mixing  length  varying 
as  e  times  tha  wall  distance;  tha  latter  imposes  a  uni¬ 
form  length  scale  at  distances  greater  than  0.13D-  from 
any  wall.  It  is  our  view  that  tha  uncertainties  "in 
initial  c<m>di cions  make  no  significant  contribution  to 
differences  between  experiment  and  computation  in  the 
hand  itself.  In  support  of  this  view,  experiments  at 
0*0*1  with  a  virtually  fully  developed  flow  at  entry  to 
the  bend  produced  ttraaarist  velocity  profiles  at  90® 
only  slightly  diffarent  from  those  shewn  in  fig.  2. 

The  eevuted  development  of  the  ft raavise  and 
radial  velocity  components  around  the  bend  is  shown 
in  figures  3  and  4.  Serious  discrepancies  between  compu¬ 
tation  and  measurement  quickly  develop.  At  45  the 
PF***e“*  radial  velocity  on  the  centreline  is  somewhat 
toe  low  and,  as  s  result,  tha  ttreswiss  calculated 
profile  at  this  position  is  biased  towards  the  inside  of 


the  bend  while  in  the  experiment  the  velocity  peak  is 
displaced  somewhat  towards  tha  outer  radius.  (It  appears 
that  tbs  experimental  measurement  of  W  ,  the  bulk  mean 
ftreamvise  velocity,  may  be  too  low  at  this  position 
causing  all  th*  measured  profiles  to  lis  above  eh*  pred¬ 
iction;  th*  differences  in  shop*  that  ar*  present  art  not 
affteted  by  this,  however).  At  90°  the  differences  ar* 
more  pronounced  including  eh*  very  strong  double  peak 
in  th*  measured  streamvis*  profile  and  its  collet* 
absence  in  th*  predictions.  The  redial  velocity 
profiles  likevis*  display  s  sharp  peak  near  th*  inner 
boundary  radius  that  is  not  reproduced  in  th*  computation. 
Similar  anomalies  are  present  at  135s  (the  experimental 
data  are  at  130°)  but  by  180°  (177°  for  measurements) 
th*  strsaafit*  profiles,  at  any  rats,  are  shoving  closer 
agreement;  the  secondary  flew  it  still  seriously  in 
error,  however. 

Clearly  something  starts  to  go  quite  seriously 
wrong  with  the  simulation  fairly  early  on  in  th*  bend  end 
a  substantially  different  flow  pattern  it  computed  at 
90s  from  that  measured.  In  search  of  the  cause  of  th* 
differences  a  number  of  adaptations  have  been  considered. 
The  most  recent  calculations  have  benefitted  from 
significant  mesh  refinement  in  th*  near-wall  region 
compared  vith  an  earlier  set  (Johnson  and  Launder  1983) 
yet  tha  differences  in  the  calculations  ar*  small  compared 
vith  the  differences  between  measurement  and  calculation. 
The  procedure  for  finding  the  secondary  wall  shear  stress 
seams  a  particular  area  of  weakness  and  so  a  test  was  made 
where  th*  wall  stresses  in  th*  x-y  plan*  were  set  to  aero. 
This  is  clearly  an  incorrect  hypothesis  but  it  served  to 
indicate  whether  wall  stress  errors  could  conceivably 
account  for  th*  large  differences.  The  largest  g 
difference  between  th*  two  practice*  occurs  at  135°  for 
which  the  primary  velocity  profiles  are  shown  in  figure 
5.  The  chiuges  produced  by  this  step  slightly  improve 
agreement  vith  experiment  and  along  2X/D„  -  0.5  a  peak 
in  V  near  the  inner  veil  is  present.  Nevertheless 
large  difference*  remain. 

It  it  well  known  that  the  k-t  model  does  not  correct¬ 
ly  capture  th*  great  sensitivity  of  real  turbulence  to 
small  amounts  of  streamline  curvature.  Although  this 
weakness  it  intrinsic  to  the  use  of  th*  Boussinesq 
stress-strain  relation,  for  tvo-dimsntionel  curved  flows 
it  has  been  found  possible  to  imitate  th*  effect  of 
curvature  on  turbulent  shear  stresses  fairly  well  by 
introducing  th*  following  term  in  place  of  th*  sink 
-cc2c2/k  in  the  c  transport  equation: 

-ct2f2  (1-0.2  Ri) 

where  Ri  *  (k/cR)zW(RH)  is  a  curvature  Richardson  nuebar 
and  R  is  th*  local  radius  of  curvature  of  a  streamline. 
This  is  given  by 

R "l  -  ((OV-VU  )*♦((*»  -VD)2*(W  -w  )k)l/R  2  (l0) 

t  t  t  t  t  t  3 /. 

where  0{  I  UU^+ViyMU^  end  *  *  (U^rV^-V2)  i 

In  places  the  secondary  velocity  field  resulting  from  this 
modification  vat  changed  by  201;  th*  effect  on  th* 
streamvis*  velocity,  however,  was  nowhere  nor*  than  3X  and 
is  thus  insignificant  compared  vith  th*  differences  her* 
in  question.  This  result  could  have  been  anticipated  for 
on*  could  not  expect  an  empirical  'fix'  on  one  stress 
compooent  in  a  two-dimensional  shear  co  be  satisfactory 
for  all  the  stress  components  in  a  copies  three- 
dimensional  flow. 

Tha  logical  next  step  in  improving  the  representation 
of  th*  Reynolds  stress  is  the  introduction  of  an  algebr¬ 
aic  stress  model  of  turbulence  (ASM)  in  place  of  the 
Boussinesq  stress-strain  relation.  Models  of  this  type 
have  been  conspicuously  successful  in  mimicking  the 
effects  of  eurvature  in  two-dimensional  shears  without 
th*  introduction  of  specially  tuned  empirical  terms 
(s.g.  Rodi  st  al,  1982).  Unfortunately,  switching  from 
a  modtl  bated  on  a  turbulent  viscosity  co  on*  where  the 
turbulent  strattet  enter  the  calculation  a*  sources  and 
sinks  is  a  severely  de-stabiliting  departure.  At  the 
time  of  writing  no  converged  results  have  been  obtained 
with  th*  QUICK  treatment  of  convection.  Vith  th*  up¬ 
wind  schema,  however,  convergence  hat  been  achieved 
though  largely,  we  believe,  because  this  approach  brings 
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ic*  on  fait t  diffusion  to  assist  stability. 

Streamwise  velocity  profiles  dram  f ran  these  results 
appear  on  figure  4  for  the  90°  position.  the 
particular  fora  of  ASM  adopted  it  that  which  reeulct 
froa  applying  todi'e  (1976)  el|tbreic  craaeport 
hypo thee ie  to  Gibaon  aad  Launder 'a  (1978)  eecoad-aoatnt 
cloaure  propoaalt.  The  Gibeos- Launder  study  which 
considered  the  caae  of  flow  past  a  tingle  plana  wall 
included  carat  rapreaencing  the  effecte  of  preteurt 
reflection  froa  the  rigid  boundary.  Here  there  are  four 
walla  praeeat  aad  chair  affect  it  attuned  to  ba  accounted 
for  by  applying  a  linear  auperpoaition;  thia  extends  to 
three— dimensions  the  usual  evo-diaeasioaal  practice. 

It  ie  clear  froa  fig.  3  chat  the  introduction  of  the 
ASM  tchaae  hat  brought  no  i^rovaaeat  in  accuracy, 
toatvhat  the  reverse.  Baaed  on  earlier  experiancee 
with  the  k-c  Bousaineaq  nodal  artificial  diffusion 
introduced  by  upwind  differencing  nodifiea  the 
streamrise  velocity  typically  by  enouats  similar  to  the 
difference  between  the  curvet  representing  ASM  aad  k-t 
predictions  in  fig.  4  (and  in  the  tent  direction) .  That 
no  i^roveatnes  are  recorded  froa  adopting  this  higher- 
level  cloaure  can  only  be  said  to  be  extraaaly  perplex¬ 
ing.  It  is  hoped  chat  by  the  cine  of  the  netting  at 
which  this  work  is  foraally  presented  an  explanation 
will  have  eaerged. 

The  heat  transfer  behaviour  considered  in  figure  6 
ie  chat  arising  from  using  che  Bousaineaq  k-t  nodal. 

The  greatest  difference  beeveen  heat  transfer  races  on 
the  inside  and  outside  of  the  bend  initially  occurs  at 
the  centre  plane.  The  secondary  flow  driving  cold 
fluid  towards  che  outer  wall  on  chit  plane  produces 
higher  beat  transfer  levels  as  is  typical  of  impinge¬ 
ment  conditions.  Conversely,  on  che  inside  of  che  band, 
che  fluid  adjacent  to  the  wall  has  arrived  chert  after 
travelling  relatively  slowly  over  a  slightly  greater 
length  (due  to  the  secondary  motion).  This  fluid  has 
thus  been  heated  considerably  more  than  chat  near  che 
oucer  wall  and  heat  transfer  coefficients  are  thus 
lower.  At  45°  the  ratio  of  che  Husaelt  nunber  on  the 
two  sides  is  3:1.  This  difference  diminishes  with 
passage  around  the  bend  as  the  secondary  aad  primary 
flow  paccarns  bacons  norm  collar;  even  at  180  however 
che  naan  level  of  Nuaselc  nuriier  is  60X  greater  on  the 
outer  than  the  inner  wall  associated  with  che  fact  that 
the  faster  moving  fluid  is  located  near  che  outer  wall 
(see  fig.  3) :  chit  produces  higher  velocity  gradients 
which  in  cum  produce  higher  kinetic  energy  level. 

Downatreaa  of  the  bend  che  secondary  flaw  dacays 
(the  mar  Inn  radial  velocity  it  only  .02  Wg  at  tea 
dianeeers)  chough  che  Huaselc  nunber  distribution 
around  the  duct  perimeter  raaaina  strongly  non-uniform. 
The  mean  level  of  Husaelt  in  che  bend  it  tone  301  higher 
than  ie  found  in  a  straight  tube  at  the  sane  Nuaselc 
amber. 

4.  CONCLUDING  BEMABKS 

Experiments  under  identical  flow  conditions  in  our 
two  institutions  using  differgnt  atthoda  of  aaaaureaant 
have  confined  that  at  the  90  poeition  in  a  square- 
sectioned  0-bend  che  etreenite  velocity  displays  a 
secondary  peak  near  the  inner  wall.  Nona  of  our  nuaari- 
cal  simulations  using  relatively  fine  meshes  aad  the 
QUICK  treatment  of  convention  haa  succeeded  in  mimick¬ 
ing  this  i^ortant  flew  feature.  Bor  the  present  we 
cannot  offer  a  convincing  exp least ion  for  this  failure: 
grid  refinement  in  the  croes-tactional  plane  haa  had 
little  influence  while  introducing  an  ASM  apparently 
brings  no  benefit.  There  are  two  seepe  we  propose  to 
take  that  any  bring  about  a  narked  improvement .  The 
first  ia  to  ceneentrata  a  far  greater  proportion  of  the 
screanwlse  calculat iooal  planes  ia  the  first  90°  of  the 
bead.  (The  present  distribution  provides  s  virtually 
uniform  spacing  around  the  bend  end  only  modest 
variations  ia  the  upstrsaa  and  downstream  tangents). 

The  second  will  be  to  abandon,  on  che  flat  walls  of  the 
bend,  the  wall-function  treatment.  Instead,  a  fine 
mash  will  be  introduced,  to  allow  coaputations  to  be 
carried  into  the  viscous  sublayer.  In  this  way  uncert¬ 
ainties  aa  to  the  appropriate  wall  boundary  condition 
are  removed.  He  have  found  froa  parallel  work  la  eha 
round-sectioned  0-bond  that  by  asauaiag  a  radial 
equilibria  pressure  distribution  across  the  viscous 


aad  buffer  regions  che  introduction  of  a  very  fine 

near -wall  grid  hat  negligible  effect  on  either  convergence 

ratsa  or  tcorage  requirements. 

With  che  velocity  field  in  such  relatively  poor 
agreement  with  data  over  note  of  the  band  no  very  defin¬ 
itive  conclusions  can  be  drawn  froa  the  detailed  heat 
transfer  behaviour.  However,  results  very  probably  give, 
in  chair  overall  pattam,  a  correct  indication  of  the 
effects  on  heat  transfer  levels:  i.e.  at  leasts  2:1  ratio 
of  heat  transfer  coefficients  on  the  outer  end  inner 
curved  walls  of  eha  bend;  a  assn  level  sons  301  higher 
chan  in  a  straight  duct  and  a  strong  non-uniformity  in 
Nusselt  number  persisting  at  least  10  disasters  down¬ 
stream. 
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INTRODUCTION 


The  movement  of  fluid  through  curved  pipes  and  bends  Is  of 
considerable  practical  and  fundamental  interest.  It  is  well 
known  that  in  such  flows  secondary  motions  arise  in  the 
cross-stream  planes  producing  a  streamwise  flow  pattern  which  may 
be  far  removed  from  that  found  in  a  straight  pipe.  The  secondary 
motions  can  be  explained  qualitatively  in  terms  of  the  response 
of  a  viscous  fluid  element  to  an  imbalance  between  the 
centripetal  acceleration  and  the  cross-stream  pressure  gradient 
Induced  by  lateral  curvature  of  the  main  flow;  see.  for  example, 
Cuming  (1957)  and  Johnston  (1978).  In  a  curved  pipe,  the  result 
is  a  cross-stream  or  secondary  motion  carrying  fluid 
symmetrically  along  the  pipe  walls  from  the  outer  to  the  inner 
line  of  symmetry  and  along  the  symmetry  plane  of  the  pipe  from 
the  inner  to  the  outer  line  of  symmetry  as  shown  in  Figure  1. 

Analysis  of  the  equations  describing  laminar  flow  through 
curved  pipes  shows  that  two  parameters  characterize  the  flow: 
the  radius  ratio  <5  =  a/R  and  the  Dean  number  De  s  6z  Re,  Berger  et 
al  (1983),  In  the  above  definitions,  a  is  the  radius  of  the  pipe 
cross-section,  R  is  the  pip£  mean  radius  of  curvature,  and  Re  is 
the  flow  Reynolds  number  W  2a/v  where  W  is  the  bulk  average 
velocity  through  the  pipe  and  v  is  the  fluid  kinematic  viscosity. 
Since  the  Dean  number  is  equal  to  the  ratio  of  the  square  root  of 
the  product  of  the  inertia  and  centrifugal  forces  to  the  viscous 
force,  it  provides  a  measure  of  the  intensity  of  the  secondary 
flow.  The  radius  ratio  6  is  a  more  direct  measure  of  the 
Influence  of  pipe  geometry  on  the  flow.  It  affects  the  balance 
of  inertia,  viscous  and  centrifugal  forces.  Berger  et  al  (1983) 
point  out  that  the  Influence  of  5  on  the  flow  through  curved 
pipes  is  not  as  well  understood  as  that  of  the  Dean  number. 

Fully  developed  flow  in  curved  pipes  has  been  the  subject  of 
extensive  research,  a  review  of  the  subject  having  recently  been 
made  by  Berger  et  (1983).  The  more  complex  entry  flow  into  a 
curved  pipe  has  been  studied  far  less  completely.  Theoretical 
attempts  to  solve  the  problem  have  been  seriously  handicapped  by 
a  lack  of  sureness  in  the  simplifying  assumptions  underlying  the 
various  analytical  approaches.  In  addition,  for  fixed  &  and  De, 
the  flow  developing  in  a  curved  pipe  is  a  function  of  the  inlet 
plane  boundary  conditions  and  these  present  some  difficulties  for 
general  values  of  Re;  see  Stewartson  et  al  (1980)  and  Berger  et 
al  (1983). 

Tao  and  Berger  (1975)  and  Stewartson  et  al  (1980)  have 
investigated  the  entry  flow  development  problem  theoretically  for 
De>>1.  In  the  former  case,  two  sets  of  equations  were  derived, 
one  for  the  invlscid  core  flow  and  the  other  for  the  three- 
dimensional  boundary  layer.  At  the  inlet  plane  zero  cross-stream 
velocity  components  and  a  uniform  axial  velocity  were  prescribed. 
Along  the  pipe  wall  zero-slip,  impermeable  conditions  were 
specified.  A  development  of  the  Karman-Pohlhausen  integral 
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method  was  used  to  solve  the  boundary  layer  equations. 

A  set  of  boundary  layer  equations,  equivalent  to  that  of 
Yao  and  Berger,  was  solved  numerically  by  Stewartson  et  al,  but 
core-f low/boundary-layer  interactions  were  neglected.  They  also 
prescribed  zero  cross-stream  velocity  components  at  the  inlet 
plane,  but  they  followed  Singh  (1979)  in  imposing  a 

potential-vortex  condition  for  the  streamwise  velocity  component. 

At  any  fixed  streamwise  location.  Re,  the  calculation  sequence 
always  advanced  from  <)>  a  0  at  the  outer  line  of  symmetry  to 
♦  a  180°  at  the  inner  line  of  symmetry.  Use  was  made  of 
preliminary  solutions  of  the  flow  field  for  all  <j>  when  R6  =  0  and 
all  R  9  when  $  s  0.  The  Keller  box  scheme  was  used  to  solve 
partial  differential  equations  for  the  leading  terms  in  a  power 
series  expansion.  The  coefficients  in  the  series  were  obtained 
seriatum  and  the  series  summed  to  yield  the  preliminary 

solutions . 

Yao  and  Berger's  solution  of  the  boundary  layer  flow  in  a 
curved  pipe  predicts  that  separation  of  the  secondary  flow  or 
circumferential  boundary  layer  will  take  place  at  a  streamwise 
location  from  the  inlet  plane  of  about  0.01a(De/5)  2  .  The  width 
of  the  separation  zone,  i.e.  the  distance  between  $  s  180°  and 
the  point  of  separation,  was  found  to  increase  with  streamwise 
location,  tending  asymptotically  to  a  maximum  value  of  59°.  This 
finding  is  in  qualitative  agreement  with  the  integral  solution 
for  fully-developed  flow  obtained  by  Barua  (1963)  who  predicts  a 
secondary  boundary  layer  separation  at  4  a  153°  ,  but  does  not 

accord  with  either  the  fully-developed  Integral  solution  of  Ito 

(1969)  or  the  numerical  studies  of  Collins  and  Dennis  (1975)  who 
do  not  predict  separation.  In  contrast,  and  for  the  first  time, 
Stewartson  et  al  predict  a  vanishing  of  the  streamwise  component 
of  skin  friction  on  the  inner  line  of  symmetry,  $  :  it  ,  at  R  e  a 
0.993  a/6*  .  The  position  of  zero  shear  stress  represents  a 
singularity  in  their  calculations  since  the  shear  stress 
increases  again  immediately  past  the  zero  point  on  the  inner  line 
of  symmetry.  From  their  study,  the  authors  concluded  that  the 
secondary  boundary  layers  "collide"  at  ^  r  180°  forming  a  radial 
Jet  which  takes  fluid  from  the  inner  to  the  outer  pipe  wall 
region  along  the  symmetry  plane. 

Experimental  measurements  obtained  by  Agrawal  et  al  (1978) 
using  the  laser-Doppler  veloclmeter  technique  in  a  transparent 
curved  pipe  with  <5  *  1/7  suggest  that  separation  of  the 

circumferential  boundary  layer,  in  the  sense  predicted  by  Yao  and 
Berger,  for  example,  may  have ,  occurred  by  R6  =  3.96  a/5*  for 
De  a  138  and  by  R9  a  6.23  a/5*  for  De  a  678.  The  values  3.96 
and  6.23  are  considerably  larger  than  those  suggested  by  either 
Yao  and  Berger's  analysis  (0.12  and  0.26)  or  the  constant  0.993 
predicted  by  Stewartson  et  al  for  the  position  of  zero  streamwise 
shear  stress.  However,  it  is  interesting  to  note  that  for  their 
higher  Dean  number  Agrawal  et  al  observed  a  striking  modification 
of  the  secondary  flow  profiles  measured  in  the  region  of  the. 
inside  of  the  bend  at  two  stations:  Re  «  1.39  a/5*  and  R6*2.31a/5*. 
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Although  the  authors  found  it  difficult  to  interpret  their 
results,  they  associated  the  phenomenon  with  some  form  of 
separation . 

In  an  effort  to  verify  the  finding  of  Stewartson  et  al, 
Talbot  and  Wong  (1982)  used  an  electrochemical  technique  to 
obtain  the  wall  shear  stress  along  the  inner  line  of  sytfcmetry  in 
a  curved  pipe  with  <5  =  1/7.  Measurements  were  made  over  the 
range  188  i  De  &  1622.  They  found  that  the  wall  shear  decreased 
to  a  minimum  value  with  increasing  De,  the  minimum  being  located 
at  R0  =0.96  a/ 5*  ,  in  close  agreement  with  the  value  predicted  by 
Stewartson  et  al.  Upstream  of  the  minimum  the  measured  wall 
shear  stress  was  also  in  agreement  with  the  predictions  by 
Stewartson  et  al,  but  on  the  downstream  side  substantial 
discrepancies  were  found  between  measurements  and  predictions. 
Berger  et  al  point  out  that  the  calculations  of  Stewartson  et.  al 
are  inaccurate  beyond  the  singularity  on  the  inner  line  of 
symmetry  and  must  be  discounted:  further  discussion  of  this  point 
is  given  by  Talbot  and  Wong  (1982). 

From  the  above  it  is  clear  that  several  important  aspects  of 
developing  curved  pipe  flow  have  not  yet  been  resolved. 
Analytical  methods  seeking  to  model  the  flow  as  an  inviseid  core 
interacting  with  a  three-dimensional  boundary-layer  lead  to 
asymptotic  results  that  are  at  variance  with  numerical 
calculations  of  the  fully-developed  form  of  the  Navier-Stokes 
equations.  Analytical  models  do  not  yet  exist  which  include  the 
effect  of  boundary-layer  separation  on  the  inviscid-core  flow 
and,  more  importantly,  the  very  existence  of  separation  is  still 
an  unresolved  issue. 

The  arrival  of  large  core  digital  computers  has  made 
feasible  the  numerical  solution  of  either  the  full  or  some 
truncated  version  of  the  Navier-Stokes  equation.  Predictions  of 
developing  laminar  flows  in  curved  pipes  by  three-dimensional 
finite-difference  procedures  have  been  reported  by  Patankar  et  al 
(1978),  Rushmore  (1975),  Liu  (1976,  1977),  Humphrey  (197?)  and 
Levy,  Briley  and  McDonald  (1083).  The  procedure  of  Patankar  et 
al  is  based  on  the  boundary-layer  equations  (with  a  uniform 
streamwise  pressure  gradient  applied  over  each  cross-sectional 
plane)  and  is  thus  applicable  only  to  pipes  with  very  large 
radius  ratios.  Levy  et  al  also  adopt  a  marching  scheme  by  taking 
the  pressure  field  from  a  potential  solution  with  a  bulk 
correction  applied  pi ane-by-pl ane  to  maintain  the  sane  mass  flow 
at  any  section.  While  this  appears  a  powerful  and  economical 
approach  to  apply  in  the  early  stages  of  development,  the  laminar 
flow  entry  problem  was  not  that  of  principal  interest  to  these 
workers  and  consequently  no  extensive  investigation  was  reported. 
Moreover,  the  approach  becomes  less  accurate  when  there  is  strong 
interaction  between  the  viscous  and  non-viscous  regions  at  high 
Dean  number.  The  other  numerical  studies  noted  above  have  been 
based  on  discretizations  of  the  full  Navier-Stokes  equations  and 
share  in  common  the  problem  of  false  diffusion  arising  from  the 
use  (for  stability)  of  upwind  differencing  with  an  inevitably 
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coarse  mesh.  As  a  result  flow  details,  such  as  small  secondary 
eddies,  tend  to  be  smeared  out;  none  of  the  studies,  for  example, 
has  reported  separation  of  the  secondary  flow. 

The  present  contribution  is  aimed  at  throwing  some  light  on 
the  various  unresolved  phenomena  discussed  above  associated  with 
developing  flows  in  pipe  bends.  Efforts  have  been  made  to  reduce 
numerical  errors  to  unimportant  levels  by  following  Pratap  and 
Spalding  (1975)  in  using  a  semi-  elliptic  rather  than  a 
fully-elliptie  treatment  (thus  permitting  considerably  finer 
meshes  than  earlier  studies),  by  adopting  the  third-order 
■quadratic  upwind  differencing  of  Leonard  (1978)  (rather  than  the 
usual  first-order  upwinding)  and  by  removing  the  numerical 
singularity  at  the  pipe  axis.  Although  numerical  resolution  is 
gradually  lost  as  the  Dean  number  is  raised,  the  results  provide 
evidence  of  secondary  flow  separation  -  albeit  very  weak  -  for  De 
greater  than  500.  They  also  show  the  occurrence  of  a  minimum 
streamwise  wall  shear  stress  on  the  inner  line  of  symmetry.  As 
the  Dean  number  is  successively  raised,  the  minimum  value 
decreases  and  occurs  at  progressively  smaller  values  of  R06*  /a 
though,  even  at  the  highest  value  of  De  considered,  this 
dimensionless  position  is  still  approximately  twice  as  far 
downstream  as  predicted  by  Stewartson's  analysis. 

During  the  documentation  stage  of  the  present  study,  the 
very  recent  thesis  of  Soh  (  1983)  came  to  our  attention.  As  in 
several  earlier  studies  he  has  used  a  fully-elliptic 
discretization.  Despite  the  resultant  coarseness  of  his 
numerical  mesh,  by  using  central  differencing  for  convective 
transport  and  by  carefully  arranging  the  mesh  non-uniformly  over 
the  domain,  his  results  reproduce  at  least  qua  1 i ta t i vel y  a 
secondary  separation.  Comparisons  are  drawn  with  Soh's  work 
wherever  possible;  the  impression  is  that  the  differences  between 
his  results  and  ours  are  attributable  to  the  finer  meshes  that 
have  been  possible  in  the  present  work. 


SUMMARY  OF  THE  NUMERICAL  SOLUTION  PROCEDURE  AND  BOUNDARY 
CONDITIONS 


The  equations  describing  the  development  of  a  viscous  fluid 
in  a  toroidal  duct  such  as  shown  in  Figure  1  may  be  written 
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In  equation  (1)  U,  V  and  W  denote  velocity  components  in  the 
circumferential  ( $  ) ,  radial  (r)  and  streamwise  directions  (9) 
respectively,  rc  =  R  ♦  aeos$  is  the  local  radius  from  the  centre 
of  the  bend,  while  in  equation  (2),  stands  for  any  of  the 
velocity  components,  the  associated  source  and  sink  terms  being 
given  in  Table  2.1.  The  operators  C(4»)  and  D(^)  are  defined  by 


C«0  =  (rrc) 


D(ip)  =  (rr  ) 
c 


■if  3 
_3$ 


3*(rcw)  +  rr(rrcv^  + 


■ill  l 
[r  3<J> 


3ip 

rcy3T 


3r 


rr 


fe(tW*>] 

"Si 


(3) 

(4) 


We  note  from  eq  (4)  that  second-derivatives  with  respect  to9 
are  omitted.  In  consonance  with  this  form  of  the  describing 
equations  a  'semi-elliptic'  finite-volume  discretization  is 
adopted,  a  description  applied  to  a  solution  where  only  the 
pressure  field  is  treated  as  elliptic.  The  pressure  thus 
requires  storing  over  the  whole  domain.  In  contrast,  the 
velocity  components  are  solved  in  a  marching  fashion  and  thus 
require  values  to  be  held  only  on  two  adjacent  r  — planes  of 
nodes.  This  type  of  solution  procedure,  first  introduced  by 
Pratap  and  Spalding  (1975),  is  particularly  attractive  in 
three-dimensional  flows  for  then  the  savings  in  memory  required 
for  the  velocity  components  (compared  with  a  fully  elliptic 
solution)  may  allow  a  sufficiently  fine  mesh  to  be  used  for  the 
pressure  field  to  reduce  numerical  errors  to  unimportant  levels. 
The  present  numerical  procedure  has  broadly  followed  the  strategy 
of  Pratap  and  Spalding,  but  several  differences  in  discretization 
and  organization  have  been  introduced  to  improve  the  numerical 
accuracy  of  the  results.  Quadratic  upstream  interpolation 
(Leonard,  1979)  is  used  to  approximate  convective  transport  in 
the  cross-sectional  plane  of  the  duct  following  the  conclusions 
of  Huang  et  al  (1983)  that  this  was  overall  the  most  accurate  of 
the  simple  treatments  of  convection  (see  also  Han  et  al,  1980). 
The  pressure/continuity  connection  is  applied  by  way  of 
Patankar's  (1980)  SIMPLER  algorithm;  in  preliminary  tests  this 
was  found  to  give  convergence  rates  an  order  of  magnitude  faster 
than  the  earlier  and  very  widely  used  pressure-correction  scheme 
SIMPLE,  Patankar  and  Spalding  (1972). 


In  the  original  semi-elliptic  scheme  of  Pratap  and  Spalding 
(1975)  no  in-plane  iterations  were  made  on  the  velocity  field  and 
thus,  of  necessity,  coefficients  of  the  difference  equations  were 
based  entirely  on  upstream  information.  Although  economical, 
this  practice  proved  to  be  inadequate  in  the  present  study  which 
has  inc lauded  tighter  bends  than  those  examined  by  Pratap  and 
Spalding  and  where,  as  a  result,  streamwise  variations  are  more 
rapid.  As  the  pressure  field,  which  was  Iterated  by  repeated 
streamwise  sweeps  over  the  solution  domain,  approached 
convergence,  first  one  and  finally  two  ln-plane  iterations  on  the 
velocity  field  were  made  with  all  coefficients  being  re-evaluated 
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using  the  current  plane  information. 

At  the  highest  Dean  numbers  a  significant  economy  in  the 
memory  required  at  any  cross-sectional  plane  was  achieved  by 
assuming  that  the  static  pressure  within  a  thin  annular  ring 
adjacent  to  the  pipe  wall  was  obtainable  from  radial  momentum 
equilibrium.  Within  this  sub-region,  which  extended  from  the 
wall  to  0.9a,  the  pressure  was  obtained  from  the  following 
degenerate  form  of  the  radial  momentum  equation: 

_1  3p  m  W2  cos$  + 

P  3r  r  r  (5) 

while  the  radial  velocity  V  was  obtained  from  the  continuity 
equation,  (4).  This  parabolic  sublayer  treatment  (PSL)  was 
originally  developed  to  facilitate  the  study  of  complex 
turbulent  flows  (Iaoovides  and  Launder,  1984)  but  it  has  also 
proved  helpful  for  the  flows  examined  here  since,  at  entry  to  the 
solution  domain,  the  boundary  layer  is  extremely  thin  and  a  fine 
near-wall  mesh  is  inevitably  required. 

There  is  no  intrinsic  need  for  a  particularly  fine  mesh  at 
the  pipe  centre.  Most  earlier  numerical  treatments,  however, 
have  put  radial  gradients  of  all  dependent  variables  to  zero  at 
r  =  0  and  to  keep  the  harmful  consequences  of  this  clearly 
incorrect  prescription  to  unimportant  levels,  a  refined  grid  has 
been  needed  in  the  vicinity  of  r  =  0.  The  central  difficulty  is 
that  at  r  *  0  (J  *  1)  there  are  many  coincident  nodes  each 

corresponding  to  a  different  circumferential  angle  4  (I).  In  the 
present  study,  to  improve  the  pattern  of  node  distribution,  the 
following  reformulation  has  been  adopted.  For  the  streamwise 
velocity  component,  W,  the  same  value  has  been  assigned  to  all 
these  coincident  centre  nodes  (W(I,1)),  this  value  being  the 
average  value  of  W  over  the  surrounding  nodes: 

NI-1 

WU.U  -  W(I,2)/ (NX-2) 

1-2 

where  1*2  and  I-  (NI-1)  correspond  with  0  and  n  radians  respectively. 
The  U  and  V  components,  however,  cannot  take  the  same  value  at 
the  centre.  Instead,  it  is  required  that  the  resultant  of  U(I,1) 
and  V( I •  1 )  should  produce  the  same  velocity  vector  irrespective 
of  the  circumferential  location.  Now  this  resultant  velocity 
must  lie  on  the  symmetry  axis:  its  value  Vres  is  obtained  as  the 
mean  of  the  radial  velocity  components  on  $  *  0  and  $  *  ir  on 
either  side  of  the  centre  nodet  The  U  and  V  velocity  components 
for  other  values  of  $  are  then  obtained  as 


-  w 
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In  fact,  Vr@8  is  the  average  of  V  on  $  -  0  and  the  negative  of  that  on  <p 


-  8  - 


ff(I,l)  -  V  sin*  ;  V(I,1)  -  V  cos* 

res  res 


At  the  tube  wall  all  three  velocity  components  are  set  to 
zero.  At  the  inlet  plane  the  streamwise  velocity  and  pressure 
are  assigned  as  uniform  over  the  plane  while  the  other  velocity 
components  are  set  to  zero.  At  the  0*180°  plane  no  constraints  are 
required  on  the  velocity  components  (which  are  treated  in  a 
boundary-layer  fashion)  while  the  streamwise  pressure  gradient 
has  been  set  uniform  over  the  section  at  a  level  needed  to 
satisfy  continuity.  This  latter  condition,  while  not  in  accord 
with  the  actual  pressure  gradient  at  180°  (which  is  affected  by 
bend-exit  effects)  has  been  found  to  affect  the  flow  pattern  only 
within  20°  of  the  exit  (for  6  =  1/7).  None  of  the  comparisons 

drawn  below  relates  to  a  position  in  the  bend  greater  than  160° 
from  the  entry. 

A  number  of  computations  reported  below  have  been  repeated 
several  times  with  different  distributions  of  nodes  in  the  three 
co-ordinate  directions  and  with  a  progressive  mesh  refinement. 
The  standard  mesh  density  employed  was  20  (radial)  x  20  (circum¬ 
ferential,*)  x  ISO  (axial,  8  )  for  the  pressure  p,  and  for  the 
velocities  28  x  20  x  2.  The  difference  in  the  number  of  radial 
nodes  for  the  velocity  and  pressure  fields  arises  from  the  fact 
that  in  the  'parabolic  sublayer'  pressure  nodes  are  not  required 
but  velocity  nodes  are. 

The  computations  have  been  made  on  a  CDC7600  computer  at  the 
University  of  Manchester  Regional  Computing  Centre.  Central 
processor  time  required  to  proceed  from  a  uniform  guessed  initial 
pressure  field  to  a  final  converged  state  where  residual  mass 
errors  summed  over  the  entire  domain  were  below  0.1*  of  the 
entering  mass  flow  ranged  from  iocibds.  for  De  »  138  to25000s.for  De  a 
2712  for  6  -  1/7. 


PRESENTATION  AND  DISCUSSION  OF  RESULTS 


Comparison  is  drawn  first  with  two  experiments  reported  by 
Agrawal  et  al  (1978),  one  in  a  bend  with  6  a  1/7  at  a  Dean  number 
of  183  (He  a  A8A)  and  the  other  with  5  a  1/20  and  De  a  565 
(Re  a  2530).  In  the  experiments,  the  bend  was  preceded  by  a 
bell-mouth  entry;  the  computer  simulation  began  with  a  uniform 
streamwise  velocity  Qand  zero  secondary  velocities  at  the  entry 
plane  to  the  bend  (0  ).  It  will  be  seen  later  from  a  comparison 
with  wall  stresses  that  the  unavoidable  mismatch  between  the 
experimental  and  computational  starting  profiles  leads  to  a  lag 
in  the  computation  by  approximately  5°  of  are  in  the  case  of 
6  *  1/7.  Streamwise  velocity  profiles  at  representative 


stations  w  i  shown  for  these  two  flows  in  Figures  2  and  3.  Due 
to  the  difference  in  refractive  index  between  the  perspex  pipe 
and  the  glycerin-water  mixture  which  provided  the  working  fluid, 
experimental  velocity  traverses  were  made  along  the  non-parallel 
lines  indicated  in  the  figures.  The  numerical  data  were 
interpolated  to  extract  velocities  along  the  same  lines. 

Near  the  bend  entry  the  potential  vortex  is  clearly  evident 
in  each  case.  For  the  lower  Reynolds  number  there  is  some 
indication  that  the  measured  boundary  layer  is  a  little  thicker 
than  the  computed,  in  part  due  to  neglect  in  the  computations  of 
any  boundary  layer  at  the  entry  plane.  As  the  shear  flow 
develops  around  the  bend,  the  secondary  motion  displaces  the 
velocity  maximum  to  the  outside  of  the  bend  and,  in  the  case  of 
the  higher  Dean  number,  producing  double  velocity  maxima  along 
some  lines.  The  computations  and  the  experiments  generally 
produce  a  strikingly  similar  behaviour.  The  somewhat  smaller 
distortion  of  the  computed  profiles  at  De  «  183  and  L/R  a  7.34 
(where  L  is  the  distance  from  the  bend  entry  measured  along  the 
circular  path  through  the  pipe  centre)  compared  with  experiment 
is,  we  believe,  probably  due  to  the  thinner  computational  Inlet 
boundary  layer.  Other  small  differences  that  a  close  examination 
reveals  are  probably  attributable  to  the  uncertainty  in  the 
measured  Reynolds  number  of  ±81!. 

The  secondary  flow  data  reported  by  Agrawal  et  al  (1978) 
were  in  fact  obtained  later  than  the  streamwise  velocities  and  at 
different  Dean  numbers,  138  and  678.  Figures  4  and  5  draw 
comparisons  between  the  measured  profiles  and  the  corresponding 
numerical  results.  In  each  ease,  the  secondary  flow  carries 
fluid  from,  the  outside  to  the  inside  of  the  bend  near  the  wall 
with  a  slow  return  flow  over  the  remainder  of  the  cross-section. 
At  the  higher  Dean  number  (or  rather,  Reynolds  number)  the 
near-wall  current  is  confined  closer  to  the  wall  due  to  the 
thinner  streamwise  boundary  layer  and  the  return  flow  pattern  is 
noticeably  more  complex  near  the  inside  of  the  bend.  *  The 
numerical  computations  mirror  the  experimental  data  reasonably 
well  but  not,  it  must  be  acknowledged,  as  well  as  in  the  case  of 
the  streamwise  profiles  discussed  above.  The  predicted  near-wall 
outer-to-inner  flow  is  thicker  than  that  measured.  This 
super f icially  might  appear  to  arise  from  numerical  diffusion  but 
the  QUICK  scheme  adopted  for  convection  is  accurate  up  to  third 
order  and  does  not  suffer  from  the  severe  numerical  smearing  to 
which  upwind  differencing  is  prone.  Moreover,  grid  refinement 
produced  negligible  changes  in  the  results  at  these  Dean  numbers. 
The  question  thus  arose  whether  the  experimental  Dean  number 
could  have  been  different  from  that  reported.  Without  claiming 
to  answer  that  question  it  is  at  least  of  interest  to  notice  that 
in  Figure  5  the  computed  secondary  flows  at  a  Dean  number  twice 
that  reported  experimentally  are  in  significantly  closer  accord 
with  the  experiment  than  the  reported  value  of  678.  At  this 
higher  Dean  number  the  secondary  profiles  reflect  the  weakening 
of  fluid  viscosity  through  the  appearance  of  secondary  maxima  and 
minima  along  a  number  of  traverse  lines  near  the  inside  of  the 


bend.  Indeed,  the  computations  at  De  s  1356  indicate  in  Figure 
5c  a  reversal  of  the  direction  of  secondary  flow  on  the  symmetry 
plane  implying  the  formation  of  a  counter  rotating  eddy  .  In  the 
experiments  the  secondary  velocity  on  the  axis  is  reduced  almost 
to  zero  along  two  traverse  lines  in  Figure  5e  but  does  not 
actually  reverse.  It  may  be  noted,  however,  that  the  measured 
secondary  velocities  appear  to  suffer  from  a  'rightward'  bias. 
That  is  to  say,  the  secondary  flow  profiles  indicate  a  mass  flow 
rate  to  the  right  that  is  from  2  to  6  times  larger  than  that  to 
the  left.  If  the  flow  were  fully  developed  in  the  axial 
direction,  continuity  would  require  precisely  the  same  rightward 
and  leftward  flow  rates  along  every  line  drawn  from  the  boundary 
to  the  symmetric  plane.  Now,  the  fact  that  the  predicted 
secondary  profiles  do  indicate,  rather  closely,  such  a  balance 
suggests  that  the  effects  of  changes  of  W  in  the  axial  direction 
do  not  make  a  major  contribution.  It  is  thus  difficult  not  to 
conclude  that  the  experiments,  for  whatever  reason,  have  given  a 
spurious  augmentation  of  the  secondary  velocity  towards  the 
outside  of  the  bend.  Thus  it  seems  probable  that  the  actual  flow 
does  Indeed  exhibit  a  secondary  flow  reversal  on  the. axis. 
Although  the  computations  -  even  those  at  De  =  1356  -  do  not 
quite  show  the  roller-coaster  appearance  of  the  experimental 
profiles,  in  view  of  the  above  discussion  the  agreement  is 
probably  satisfactory.  Agrawal  et  al  (1978)  speculated  that 
'separation'  of  the  secondary  flow  may  have  occurred  at  their 
data  collection  point  closest  to  the  inside  of  the  bend.  In 
fact,  though  one  cannot  distinguish  it  in  the  figure,  the 
computed  wall-ad jacent  velocities  at  this  position  do  take  very 
small  negative  values  as  may  be  inferred  from  the  variation  of 
circumferential  shear  stress  in  Figure  12. 

A  clearer  impression  of  the  secondary  flow  pattern  is 
perhaps  conveyed  in  Figures  6  and  7  which  show  secondary  flow 
'vectors'  at  different  axial  stations  for  two  values  of  De.  At  a 
Dean  number  of  138  the  ■pattern  of  the  secondary  flow  streamlines 
changes  little  as  the  flow  develops  around  the  bend  but  the 
magnitude  of  the  secondary  velocities  in  different  regions 
changes  considerably,  the  maximum  levels  being  reached  at  about  6 
radii  into  the  bend  (50°).  At  De  *  678  the  initial  simple 
cellular  pattern  evolves  into  a  kidney  shape  by  12.8  radii 
downstream  with  very  rapid  spatial  variations  in  magnitude  and 
direction  of  the  secondary  flow.  In  marked  contrast  with  the 
situation  in  Figure  7b,  the  secondary  velocities  on  the 
centreplane  are  nearly  zero.  Soh  (1983)  reports  a  qualitatively 
similar  behavlous  at  this  Dean  number.  The  behaviour  for  a  Dean 
number  of  1356  (not  shown)  is  similar  to  that  at  678  except  that 
the  principal  return  flow  is  squeezed  somewhat  closer  to  the  pipe 


^  interesting  to  note  that  Azzole  end  Humphrey  (1984)  have  measured 
such  a  counter  rotating  eddy  near  the  symmetry  plane  for  turbulent  flow 
in  a  180°  bend 


boundary  leaving  a  larger  core  with  virtually  no  secondary 
aotion. 

It  is  of  interest  to  return  briefly  to  the  streamwlse 
velocity  component  to  observe  how'  its  isovels  are  distorted  by 
the  secondary  flow  field.  Figure  R  presen tsQ measured  and 
predicted  contours  for  a  Dean  number  of  565  at  83  from  the  bend 
entry.  Soh  (1983)  provides  computed  results  for  this  flow 
condition  at  the  same  location  and  his  predictions  are  included 
in  the  figure.  The  contour  plots  convey  a  clear  Impression  of 
how  the  secondary  flow  (similar  to  that  shown  in  Figure  6c)  pulls 
out  the  axial  contours  as  fluid  flows  along  the  walls  from  the 
outside  to  the  inside  of  the  bend.  Because  the  bulk  of  the 
return  flow  also  takes  place  along  the  periphery  of  the  pipe, 
however,  the  profiles  are  folded  back  on  themselves  forming 
hook-like  contours  or "fingers".  The  present  numerical  results 
mimic  closely  the  measurements  save  that  the  computed  fingers  are 
a  little  wider.  Soh • s  coarse  grid  calculations  yield  contours 
for  W/W  equal  to  1.2  and  1.4  enclosing  smaller  regions  of  the 
flow  than  either  the  measurements  or  the  present  calculations; 
the  main  features  of  the  flow  are  nevertheless  quite  well 
predicted . 


Figures  9-12  relate  to  the  distributions  of  wall  shear 
stress  around  the  inside  of  the  bend,  a  topic  that  has  been  the 
main  concern  in  the  papers  by  Stewartson  et  al  (1980)  and  Talbot 
and  Wong  (  1982).  Figure  9  shows  the  development  of  the  axial 
wall  shear  stress  along  the  inner  line  of  symmetry  at  four  Dean 
numbers;  also  included  is  the  behaviour  predicted  by  Stewartson 
et  al.  As  discussed  in  the  Introduction,  the  analysis  developed 
by  these  workers  gives  a  vanishing  shear  stress  at  Z  =  6(R/a)*  - 
0.943  which  represents  a  point  of  singularity  in  the  solution 
since  immediately  downstream  therefrom  the  shear  stress  rises 
sharply  then  approaches  monotonlcally  an  asymptotic  value. 
Stewartson  et  al  comment  that  their  analysis  is  strictly 
applicable  only  for  very  small  values  of  S  and  for  De  >>  1. 

Certainly,  as  the  Dean  number  is  successively  raised  the  present 
numerical  results  shift  in  the  direction  of  that  limit,  i.e.  the 
minimum  dimensionless  wall  shear  rate  falls  as  De  Increases  and 
the  minimum  value  occurs  at  progressively  smaller  values  of  Z 
(though  even  at  De  *  2712  the  minimum  is  reached  about  twice  as 
far  downstream  as  the  predicted  singular  point).  Downstream  of 
the  minimum,  the  numerical  solutions  display  a  damped  oscillatory 
behaviour,  the  amplitude  growing  as  the  Dean  number  is  raised. 
This  behaviour  is  evidently  in  striking  contrast  with 
Stewartson's  result,  yet  is  at  least  in  qualitative  agreement' 
with  the  invlscld  analysis  of  Hawthorne  (1951). 

It  would  have  been  Interesting  to  extend  the  numerical 
results  to  higher  Dean  numbers  but  this  was  not  feasible  since  to 
achieve  sensible  grid  independence  for  larger  De  than  those 
reported  would  have  required  finer  meshes  -  and  thus  more 
in-core  storage  than  was  available  to  us.  It  is  of  interest  to 
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note  the  effect  of  grid  refinement  on  the  present  solutions. 
Figure  10  shows  results  obtained  for  the  three  highest  Dean 
numbers  with  the  standard  20  x  28  x  150  grid  and  with  a  coarser 
version:  20  x  20  x  100.  At  Dean  numbers  of  678  and  1356  the 
changes  in  shear  stress  arising  from  grid  refinement  are  rather 
small  and  the  .trend  is  generally  to  raise  the  shear  stress 
slightly.  The  change  is  more  substantial  at  the  highest  Dean 
number,  the  minimum  shear  stress  being  raised  by  a  factor  of  3. 
It  is  ironic  that  for  this  last  case  the  coarser  grid  gives  a 
behaviour  up  to  the  position  of  minimum  shear  stress  much  closer 
to  that  predicted  by  the  analysis  of  Stewartson  et  al. 
Nevertheless,  the  finer  grid  solution  brings  the  behaviour  at 
this  highest  Dean  number  much  more  into  line  with  the  numerical 
results  at  lower  values  of  De. 

The  conclusion  that  the  solution  of  Stewartson  et  al  did  not 
adequately  describe  the  flow  downstream  of  the  singularity  was 
originally  drawn  by  Talbot  and  Wong  (1982)  on  the  basis  of 
experimental  shear  stress  data  obtained  by  an  electrolytic 
method.  Comparisons  with  these  measurements  are  shown  in  Figure 
11  where  the  present  computational  curves  are  all  displaced  to 
the  left  a  dimensionless  distance  0.2  (corresponding  to  9.3°  of 
arc),  an  arbitrary  adjustment  to  try  and  account  for  the  effects 
of  the  inlet  contraction.  For  De  s  678  the  computed  curve 
corresponding  to  the  streamwise  wall  stress  at  8ir/9  is  also 
included  to  allow  comparison  with  the .  resul tant  stress  along  this 
line  measured  by  Choi  et  al  (1979)  .  Agreement  between  the 

experimental  and  numerical  results  is  somewhat  mixed.  At  a  Dean 
number  of  183  the  computed  values  are  some  20%  below  the  data, 
while  at  the  higher  Dean  numbers  the  cluster  of  data  points 
around  Z  =  1.0  give  substantially  lower  values  than  predicted. 
It  is  hard  to  ascribe  a  level  of  accuracy  to  the  experiments: 
the  calibration  curve  from  Talbot  and  Wong  suggests  that  the 
stress  levels  are  systematically  low  by  an  amount  ranging  from 
15-25%  depending  on  the  surface  strain  rate,  though  no  estimates 
of  other  uncertainties  are  provided.  Apart  from  the  case  of  the 
lowest  Dean  number  the  impression  conveyed  by  the  data  seems  to 
be  that  they  scatter  about  the  numerical  predictions  rather  than 
display  conclusive  differences.  Talbot  and  Wong  inferred  from  a 
comparison  of  their  measurement  with  those  at  8*/9  from  Choi  et 
al  (1979)  that  the  circumferential  wall  shear  stress  at  this 
position  was  much  smaller  than  the  streamwise  stress 
conclusion  which  conflicted  with  the  predicted  behaviour  o 
Stewartson  et  al.  The  present  study  provides  strong  support  for 
Talbot  and  Wong’s  conclusion.  The  circumferential  stress  along 
8x/9,  shown  in  Figure  12,  is  an  order  of  magnitude  smaller  than 
the  streamwise  component  except  in  the  vicinity  of  its  maximum 
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Because  the  circumferential  stress  at  this  position  is  small  compared  with 
the  axial  stress,  the  resultant  stress  does  not  differ  from  the  axial  value 
by  more  than  12 
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value.  We  note  that  the  streamwise  variation  of  this  component 
is  essentially  Independent  of  Dean  number  as  far  as  Z  s  1.5. 
Moreover,  weakly  negative  values  of  circumf erential  shear  stress 
occur  for  De  s  678  in  the  range  2.74Z49.8  and  for  De  =  1356  in 

the  range  3.7-SZS5.9.  At  the  highest  Dean  number  the  secondary 
shear  stress  remains  positive,  though  close  to  zero  for  Z>  U.  It 
ought  to  be  said  that  this  last  result  is  not  conclusively 
established  since  it  is  possible  that  a  further  major  grid 
refinement,  while  producing  negligible  changes  to  the  streamwise 
or  secondary  velocity  field,  could  nevertheless  change  the 
circumferential  stress  from  a  very  weak  negative  value  to  an 
equally  weak  positive  value  -  or  vice  versa. 


CONCLUSIONS 


Careful  numerical  solutions  have  been  obtained  of  several 
laminar  flows  developing  in  180°  bends  of  circular  cross-section 
that  have  been  the  subject  of  laser  Doppler  studies  by  Talbot  and 
his  colleagues.  Given  the  small  but  unquantifiable  mismatch 
between  the  computational  and  the  experimental  entry  conditions, 
agreement  between  computed  and  measured  streamwise  profiles  for 
Dean  numbers  of  183  and  565  are  in  extremely  close  agreement. 
There  is  less  complete  agreement  with  the  secondary  velocity 
profiles  and  one  possibility,  that  would  he  consistent  with  the 
present  results,  would  be  that  the  experimental  data  (which  were 
obtained  in  a  separate  study  from  the  streamwise  profiles)  were 
at  a  higher  Dean  number  than  reported. 

The  present  computational  results  indicate  a  gradual 
approach  towards  the  initial-region  behaviour  predicted  by 
Stewartson  et  al  (1980)  as  the  Dean  number  is  raised.  Even  at  a 
Dean  number  of  2712,  however,  there  are  still  marked  differences 
from  Stewartson's  solution.  Downstream  of  the  point  of  minimum 
shear  stress  on  the  inside  wall  the  numerical  results  indicate  an 
oscillatory  development  of  the  streamwise  wall  stress,  the 
overshoot  increasing  as  the  Dean  number  is  raised.  This 
behaviour,  which  is  at  least  in  qualitative  agreement  with  the 
data  of  Talbot  and  Wong  (1982),  and  inviscid  flow  calculations  of 
the  secondary  flow  by  Hawthorne  (1951),  is  in  striking  contrast 
with  the  predictions  of  Stewartson  et  8l  (1980)  which  show  a 
monotonic  approach  to  steady  state  conditions. 
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Figure  5  Secondary  velocity  profiles 
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•  The  Computation  of  Momentum  and  Heat  Transport  in 
j  Turbulent  Flow  around  Pipe  Bends 

j  H.  Iacovides  and  B.E.  Launder 


A  nvmerical  solving  procedure  is  described  based  on  a  semi- 
elliptic  discretization  of  the  averaged  equations  of  motion 
describing  turbulent  flow  through  curved  circular- sectioned 
tubes.  It  adopts  Leonard's  QUICK  treatment  of  convection 
and  Patankar's  SIMPLER  algorithm  for  handling  the  pressure- 
velocity  connection.  Applications  are  reported  of  laminar 
and  turbulent  flew  in  a  90°  bend  (adopting,  in  the  latter 
case,  the  standard  k-e  Boussinesq  viscosity  model  for  the 
Reynolds  stresses) .  Solutions  of  the  thermal  energy 
equation  for  this  case  indicate  a  marked  rise  in  the  averane 
heat  transfer  coefficient  at  any  section  but  with  a  five¬ 
fold  variation  in  Nu  between  the  inside  and  outside  of  the 
tube. 


1 .  Introduction  | 

» 

Turbulent  flow  around  pipe  bends  is  a  common  feature  of  heat 
exchangers,  condensers,  boilers  and  other  heat-transfer  equipment. 
iSuch  flows  also  arise  in  the  internal  cooling  passages  of  various 
components  of  power  generation  equipment  from  gas  turbine  blades 
to  the  rotating  cores  of  electrical  generators.  Beside?  the  very 
complex  three-dimensional  flow  that  occurs  in  the  bend  itself,  the 

secondary  flow  generated  in  the  bend  can  carry  over  to  have  an 

appreciable  effect  on  the  velocity  and  temperature  fields  many 
diameters  downstream.  ! 

! 

Despite  the  importance  of  the  flow,  little  is  known  of  the 
detailed  structure  in  the  turbulent  flow  regime.  For  the  most 
;part,  studies  have  reported  overall  pressure-loss  coefficients  for 
ibends  of  different  angles  and  different  ratios  of  pipe  diameter  to 

bend  radius.  Exceptions  are  the  measurements  of  mean  velocity  by 

Rowe  (1970)  in  a  rather  gentle  180°  bend  and  the  recent  experiment 
lof  Enayet,  Gibson,  Taylor  and  yianneskis  (3982)  which  reports  the 
^development  of  the  mean  and  rms  turbulent  streamwise  velocities 
jaround  a  'tight'  90°  bend  (pipe  diameter:  mean  radius  of  curvature 
|»  2.8:1).  There  are  no  local  heat  transfer  data  available,  so  far 
•as  we  know,  in  the  development  region  of  a  circular  sectioned 
ipipe,  though  a  joint  programme  of  experiments  on  both  the  flow  and 
thermal  fields  is  now  in  progress  at  IJMIST  and  the  University  of 
California,  Berkeley  (Professor  J.A.C.  Humphrey).  j 

i 

The  advance  in  speed  and  core  capacity  of  computers,  coupled 
with  improvements  in  numerical  methods  for  solving  the  fluid  flow 
equations,  make  it  feasible  to  undertake  computer  simulations  of 
.three-dimensional  flows;  these  can  provide  far  more  detail  of  the 
^velocity  and  thermal  fields  than  could  any  experimental 
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realization  of  “the- "same-’ flow7  Whether  or'  not  the  resulting 
calculations  provide  a  close  approximation  of  the  real  flow 
depends  on  how  small  the  com pu tor  has  been  able  to  keep  his. 
numerical  error  and  on  how  well  the  unknown  turbulent  momentum  andf 
heat  fluxes  are  represented  by  the  'turbulence  model'  adopted.  In 
the  case  of  developing  turbulent  flow  around  pipe  bends,  the  first, 
computations  were  reported  by  Pratap  and  Spalding  (1975).  The 
paper  introduced  a  novel  simplification-  of  the  full 
three-dimensional  Navier  Stokes  equations  in  that,  while  the  fully 
elliptic  form  of  the  pressure  equation  was  retained,  a  parabolic 
'truncation  was  adopted  for  the  velocity  field  wherein  streamwise 
momentum  diffusion  was  omitted.  This  simplification  allowed  the 
velocity  field  to  be  stored  on  only  two  adjacent  streamwise  planes 
instead  of  over  the  whole  domain.  The  savings  in  core  that  such 
semi-elliptic  schemes  allow  enables  a  far  finer  mesh  to  be 
employed  for  the  pressure  field  (which  must  be  stored  over  the 
full  domain)  than  would  otherwise  be  possible.  Pratap  and 
Spalding  (1975)  achieved  moderately  good  agreement  in  their 
simulation  of  Rowe's  (1970)  experiment,  though  some  of  the 
ifeatures  of  their  original  calculation  method  are  unsuitable  for 
computing  the  practically  far  more  important  case  of  flow  around 
'sharp  bends  with  R/d  in  the  range  2-4.  : 

J 

Perhaps  the  most  serious  limitation  of  the  Pratap-Spalding 
scheme  for  such  flows  is  their  use  of  "wall  functions"  (c.f. 
;Launder  and  Spalding,  1974),  based  on  experimental  data  for  plane 
jtwo-dimensional  turbulent  flows,  to  provide  the  near-wall  boundary 
conditions  on  the  velocity  field.  This  approach  provides  an 
'unreliable  basis  for  fixing  the  secondary  velocity  parallel  to  the 
Wall  and  it  is  difficult  to  conclude  other  than  that  the  only 
satisfactory  treatment  is  to  carry  the  integration  to  the  wall. 
(This  approach  has  been  adopted  by  McDonald's  group  (McDonald, 


1982)  in  computing  turbulent  flow  around 


bend  of  square 


cross  section  as  part  of  their  contribution  to  the  1981  Stanford 
Conference  on  Complex  Turbulent  Flows  (Kline  et  al,  19R2)  .  The 
results  of  this  computation,  while,  overall,  not  superior  to  those 
of  other  submissions,  did  nevertheless  give  a  more  accurate 


of  other  submissions,  did  nevertheless  give  a  more  accurate 
(account  of  the  development  in  the  first  30°  of  the  bend  where 
,there  is  a  strong  secondary  flow  confined  very  close  to  the  wall 
Carrying  fluid  from  the  outside  to  the  inside  of  the  bend.  The 
jnumerical  procedure  developed  by  McDonald's  group  obtains  the 
pressure  distribution  as  a  one-dimensional  correction  to  a 
potential  flow.  This  simple  approach  works  very  well  when  the 
boundary  layers  are  thin  but^appears  to  be  a  major  factor  in  the 
deterioration  in  their  predictions  of  the  square  bend  experiment 
beyond  45°.  A  further  potential  source  of  error  is  the  use  of 
the  very  simple  mixing-length  hypothesis  to  compute  the  Reynolds 
{stresses.  The  scheme  entirely  neglects  transport  effects  on 
turbulence;  moreover,  as  Levy  et  al  (1983)  comment,  the 
prescription  of  a  boundary-layer  thickness  -  required  to  fix  the 
mixing  length  -  becomes  arbitrary  in  such  a  three-dimensional 
flow. 


I  The  present  contribution  reports  our  attempts  at  achieving  a 
jnore  satisfactory  numerical  and  physical  computational  model  of 
turbulent  flow  around  bends  than  has  hitherto  been  reported.  We 
adopt,  following  Pratap  and  Spalding  (1975),  a  semi -el  1  i pt i c 
{approach  though,  like  Levy  et  al  (1983),  we  eschew  the  use  of  wall 


i  functions  and  instead  extend  the  computations  to  the  pipe  wall 
using,  in  the  low-Reynolds-number  sub-layer,  the  Van  Driest  (1956) 
version  of  the  mixing  length  hypothesis.  In  the  main  flow, 
however,  the  k~e  Boussinesq  viscosity  model  is  adopted  to  allow 
some  account  to  be  taken  of  transport  effects  on  turbulence.  i 

I 

Extensive  comparisons  are  drawn  with  the  laminar  and  turbulent 
flow  data  of  Enayet  et  al  (19B2)  in  a  90°  bend.  We  report  also 
the  computed  behaviour  for  the  thermal  field  for  this  geometry. 
Finally,  comparisons  are  drawn  with  the  Seban  and  McLaughlin 
(1963)  measurements  of  local  Nusselt  number  in  fully  developed 
flow  through  coiled  tubes.  i 

J 

2.  Summary  of  the  Computational  Procedure  i 


2.1  The  Describing  Equations  and  Boundary  Conditions 


The  flow  is  analysed  in  the  toroidal,  r,  0,  e  system  shown  in 
figure  1.  If  we  adopt  the  Boussinesq  turbulent  viscosity 
concept  to  represent  the  momentum  and  heat  fluxes  due  to  the 
turbulent  motion,  the  describing  continuity,  momentum  and' 
energy  equations  may  be  written  as  follows:-  i 


Continuity 


!jxcU  +  l7rrcV  +  fe*tW  0 


iMomentum  and  Enerc 


P(C(*)  +  Sc(i|0)  -  D0J0  +  Sd(i|.)  +  spW 


(2.1) 


(2.2) 


In  equation  (2.1)  U,  V  and  W  denote  the  velocity  components  in 
the  0,  rr6 directions,  rc is  the  local  radius  from  the  centre  of  the 
bend  (R+rcosG  ),  while  in  equation  (2.2)  V  stands  for  any  of  the 
velocity  components  or  the  temperature.  The  associated  source  and 
sink  terms  are  given  in  Table  2.1  and  the  operators  C(  <i>  )  and  D 
(♦)  denote: 


«♦>  2  <rrc)-‘  yj- 


rcU«0  +  ^-(rrcV<|<)  + 


“ 

h (rwW. 


•D<*>  5  (rrc>[H?(tc(!r*  ^t)lf'). 


*  I 

The  quantity  ov  is  unity  in  all  but  the  energy  equation  where  it 


v  io  uni  Ljr  ah  an  mul  uie  cncty  y  c^ua  liuii  wiicic  it 

denotes  the  molecular  Prandtl  number  of  the  fluid.  Notice  that 
second  derivatives  with  respect  to  e  are  dropped,  an  essential 
feature  of  semi-elliptic  treatments.  Where  computations  are  to  -be 
'made  in  straight  approach  sections  of  the  pipe  or  in  downstream 
tangents,  the  mean  radius  of  curvature  R  is  set  to  some  suitably 
llarge  value  in  these  regions. 


Sc<*> 


sp«0 


V*> 


VU  4  tW2sin»} 
rc 


I  lE 

7  30 


7cfe^Veff^^v»*77>Veff(!T-u)' 
+  ueff  {r(7}^«S-  ^750-^ffVf  )J 

2u  f-*in0  3W 

-  —  ■  -. .  {Usin«  -  Vcost  -  —  } 
r  2  30 


cos0W2 


-  {• 


-«?> 


-  { 


♦  {- 


sin«UV 

c 

CO«0VW 


lE 

3r 


I  iE 

r  36 
c 


\eemvt  *  " 


13,  3  ,0.  .  1  3  ,  3V 

7  30  *  eueff3r  (7*  }  7r  3r  trrcuef£  3r 
c  c 


-*"^££  (—  ♦  v>  +  A  —  Irp  —  <H» 

r  30  r  30  *  3r  r 

c 

2v  ffcos0  3W 

♦  —Si i -  (UsinO  ♦  VcosO - } 

r  30 


7r  k  {u*ff(lr +  Wsin?)} 

c 

*  7rc  fc  {2Wl?  '  Wcos,)} 

♦-2  h  {u.ff(H" 2Usin9  *2Vcos*)} 

rc 

u  „sin0  3  W  3  W 

- 7“  3?  (7>  *  ueffco5*7r  (7.) 

c  c 

V#f£  r  tin  0  3  U  u^ffcos0  3V 

T  30  (rJ  “T  (33) 

C  C 


Table  2,1  Source  and  diffusion  terms  in  mean  flow  equations 
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i  ~'Calculations'~are  "made  over  the  semi-circular  half  cross; 
I  section  bounded  by  the  diametral  plane  of  symmetry  passing  through 
the  centre  of  the  pipe  bend.  The  velocity  normal  to  that  plane  is 
set  to  zero,  while  the  gradients  of  W  and  V  normal  to  the  plane  are 
zero.  Along  the  pipe  wall  no-slip  conditions  are  applied.  The. 
velocity  components  and  pressure  are  all  prescribed  at  entry  to 
the  flow  domain  while  at  exit  the  second  derivative  of  the  static 
pressure  in  the  streamwise  direction  is  set  to  zero;  the  velocity 
and  energy  equations  are  only  first  order  in  6  and  no  outlet 
boundary  condition  is  required.  j 

As  described  in  detail  elsewhere  (Launder,  Johnson  and 
Iacovides,  1983;  Humphrey,  Iacovides  and  Launder,  1984)  a  new 
treatment  of  the  U  and  V  boundary  condition  at  the  tube  axis  has 
been  evolved  in  this  work  that  removes  the  usual  singularity 
problem.  The  need  for  special  care  arises  because  the  axis 
coincides  with  the  location  of  many  U  and  V  nodes  each  corresponding 
to  different  values  of  the  co-ordinate  0  .  Now,  since  the  nodes  in 
question  lie  at  the  symmetry  plane,  we  know  that  for  every  such 
node  the  resultant  velocity  in  the  plane  of  the  cross  section  must 
be  the  same  and  directed  along  the. symmetry  plane  -  since  as  noted 
in  §2.1.  there  can  be  no  velocity  component  normal  to  this  plane. 
We  thus  conclude  that: 

UCI,1)  -  Vrefl  sin0  (I) 

V(I»1)  -  V  cos  0(1) 

res  (2.5) 

where  Vreg  is  the  average  radially  directed  velocity  at  the  axis,' 

K'Vo-  V.>- 

Application  of  these  conditions  has  removed  the  need  for  a 
very  small  near-axis  cell  adopted  in  earlier  treatments  (e.g.  Levy 
et  al,  1983)  . 

i 

I 

2.2  The  Turbulence  Model  j 

1 

Although  in  the  future  the  authors  aim  to  include  a  model  of 
turbulence  based  on  closure  of  the  stress  and  heat-flux  transport 
equations,  the  present  computations  all  adopt  the  Boussinesq 
concept  of  an  isotropic  turbulent  viscosity  y  ;  in  the  momentum 
equations  is  unity,  while  in  the  energy  equation  it  represents 
the  turbulent  Prandtl  number  and  is  assigned  the  uniform  value 
0.9.  j 

It  cannot  be  said  in  advance  how  serious  the  undoubted 
conceptual  weakness  of  such  an  approach  will  in  fact  prove  to  be. 
It  seemed  desirable  to  us  to  adopt  this  .  simple  hypothesis 
initially  since  its  use  facilitates  the  numerical  solution  of  the 
describing  equations;  moreoever,  it  will  be  desirable  to  have  a 
'baseline'  prediction  available  against  which  to  judge 
computations  to  be  obtained  later  with  more  elaborate  models  of 
the  turbulence  field.  j 

Prom  the  wall  extending  to  a  radius  of  0.9*a  -  0. 9*ai,where  a 
is  the  pipe  radius,  the  turbulent  viscosity  is  given  by  the* mixing 
length  hypothesis  (Prandtl,  19?*}  extended  to  three-dimensional 
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flows  by' requiring  "consistency  with  the  local  equilibrium  form  of 
the  turbulent  kinetic  energy  equation.  In  the  near-wall  zone/ 
gradients  with  respect  to  r  far  outweigh  other  terms  and  so: 


i 


(2.6) 


that 


The  mixing  length  lm  is  given  by  Van  Driest's  (1956)  proposal 


1  -  <y(l  -  exp  -  j/Ut/26v) 


(2.7) 


where  y  is  the  distance  from  the  wall  (a -r)  and  UT  is  the  friction 
velocity  based  on  the  resultant  wall  shear  stress.  The  von  Karman 
constant  <  takes  the  value  of  0.419.  j 

Over  the  main  region  of  flow  the  turbulent  viscosity  is' 
obtained  from  the  standard  high  Reynolds  number  version  of  the  e 
Boussinesq  viscosity  model  (Jones  and  Launder,  1972,  Launder  and 
Spalding,  1974)  wherein  i 

y  -  c  px2/e 

t  u  (2.8) 

|and  the  turbulent  kinetic  energy  k  and  its  rate  of  dissipation 
are  found  from  the  solution  of  transport  equations  of  the  same 
type  as  equation  (2.2)  in  which  the  source/sink  terms  take  the 
following  form: 
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The  empirical  coefficients  appearing  in  the  k  and  c 
'equations  take  the  following  standard  values: 


C  -  .09;  C  -  1.44;  C 

VI  2 


1.90;  -  1.0;  ae  -  1.22 
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At  the  interface  with  the  mixing  length  hypothesis,  the  values' 
of  k'  and  e  are  fixed  by  requiring  that  the  viscosity  given  by 
equations  (2.6)  and  (2.8)  should  be  the  same  and  by  taking  the 
turbulent  length  scale  l(=k 3/2  / e)  equal  to  l  / cy3/4  as  in  a  simply 
sheared,  equilibrium,  near-wall  flow.  Theresultant  expressions 
are 

k  -  u*  /c  y*  ;  e  -  c  k2/p 

t  u  m  y  c 

where  M t  is  the  value  of  turbulent  viscosity  given  by  eq  (2.6). 

0 

2.3  Discretization 

| 

The  present  computational  procedure  has  taken  as  its  starting 
point  a  code  for  solving  the  f ully-ell iptic  laminar  flow  equations 
in  toroidal  co-ordinates  kindly  provided  by  Professor  J.A.C. 
Humphrey  (Humphrey,  1977)  .  This  had  been  developed  within  the 
framework  of  the  well-known  TEACH  family  of  computer  programs. 
Like  other  finite-volume,  primitive  variable  procedures,  the 
differential  equations  presented  in  §2.1  and  §2.7  are  discretized 
by  integrating  them  over  small  contiguous  control  volumes  which 
together  cover  the  whole  flow  domain.  To  each  control  value  is 
attached  a  discrete  value  of  the  dependent  variable.  A  staggered 
arrangement  of  grid  nodes  is  adopted  with  velocity  components 
located  at  the  boundaries  of  the  control  volumes  for  the  scalar 
variables  (p,T,k  and  e  ).  We  retain  the  usual  TEACH  practice  of 
assuming  a  linear  variation  of  dependent  variable  between  nodes  in 
evaluating  diffusion  processes  while  treating  source  and  sink 
terms  as  uniform  over  each  control  volume  with  a  value  equal  to 
that  at  the  node. 

In  the  treatment  of  convection,  the  upwind/central  hybrid 
differencing  normally  incorporated  in  TEACH  is  retained  only  for 
the  streamwise  velocity  component  W.  For  the  components  in  the 
plane  of  the  cross  section,  the  more  accurate  quadratic  upwind 
interpolation  (Leonard,  1979)  has  been  incorporated.  (See  Han  et 
al,  1981  and  Huang  et  al,  1983  for  a  comparison  of  the  performance 
of  quadratic  and  hybrid  differencing) . 

i 

Pratap  and  Spalding  (1975)  adopted  the  SIMPLE  procedure 
(Patankar  and  Spalding,  1972)  for  use  in  their  semi-elliptic 
computations.  We  have  found,  however,  that  the  newer  scheme 
(SIMPLER  (Patankar,  1980)  produced  rates  of  convergence  an  order  of 
magnitude  faster.  The  original  pressure-velocity  iteration 
sequence  has  had  to  be  modified  somewhat  to  allow  it  to  fit 
compatibly  within  the  present  semi-elliptic  treatment  with  only 
the  pressure  available  in  three-dimensional  storage:  (i)  The 
solution  begins  with  a  guessed  initial  pressure  field  starting 
from  the  initial  (upstream)  plane  and  marching  plane  by  plane 
downstream.  At  each  step  with  the  pressure  field  fixed,  the  U  and 
V  momentum  equations  are  first  solved  on  a  plane  followed  by  that 
for  the  W  velocity  (located  one  half  cell  downstream).  (ii)  The 
velocity  field  is  next  corrected  by  way  of  the  standard 
pressure-correction  equation  in  order  to  satisfy  continuity. 
j(iii)  The  current  plane  pressure  field  is  then  updated  by  solving 
ithe  Poisson  equation  for  pressure _oyer_.the_pl.ane  using  the  updated 
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'velocities.'  '  "  (iv)  The  currerTt  plane  velocity  field  is 
re-calculated  using  the  updated  pressure  field.  (v)  The  current 
plane  velocities  and  the  pressure  field  on  all  downstream  planes 
are  adjusted  by  way  of  the  Pratap-Spalding  bulk  pressure 
correction  to  ensure  that  the  prescribed  mass  flow  rate  is  passing 
that  section.  (vi)  The  current  plane  velocity  field  is 
re-corrected  via  the  pressure-correction  equation.  (vii)  For 
turbulent  flow,  the  turbulent  kinetic  energy  and  dissipation  rate 
are  solved  and  the  turbulent  viscosity  is  updated. 

The  above  steps  are  repeated  at  each  succeeding  plane  and, 
when  a  complete  downstream  sweep  has  been  completed,  the  sequence 
begins  over  again  at  the  upstream  plane  unless  the  solution  has 
fully  converged. 

As  the  calculation  approaches  convergence,  first  one  and 
later  two  iterations  are  introduced  at  each  plane  in  solving  the 
momentum  equations  in  order  to  re-form  coefficients  and  source 
terms  in  terms  of  current  plane  values.  (When  marching  without 
iteration,  one  necessarily  has  to  evaluate  these  quantities  from 
upstream  values).  Although,  apparently,  Pratap  and  Spalding 
(1975)  did  not  incorporate  such  iterations  into  their 
computations,  the  bend  flows  considered  in  their  study  were  far 
milder  than  those  of  the  present  study  and  thus  streamwise 
variations  were  less  rapid. . . 

Within  a  thin  annular  ring  extending  from  the  wall  to  o.P5a 
a  much  simpler  and  economical  numerical  procedure  has  been  adopted 
in  our  most  recent  calculations.  The  region  is  treated  as  a 
parabolic  sublayer  (PSL) ,  lacovides  and  Launder  (19P4)  in  which, 
[for  a  given  0  and  0  ,  the  pressure  at  any  point  in  the  sublayer 
above  that  at  the  first  node  outside  the  sublayer  is  obtained  by 
assuming  radial  equilibrium.  Thus,  no  pressure  nodes  are  required 
in  the  PSL,  a  feature  that  allows  the  use  of  a  fine  near-wall  mesh 
for  the  velocity  field  without  prohibitive  core  demands  (the 
pressure,  it  will  be  recalled,  is  the  only  variable  for  which 
three-dimensional  storage  is  needed) .  A  further  simplification  is 
that  the  radial  momentum  equation  is  not  solved,  the  V  velocities 
being  obtained  directly  by  application  of  continuity  to  the 
control  volume  surrounding  the  pressure  node.  j 


3.  Numerical  Results 


UMRCC'S 

purely 

There 

three- 


The  computations  reported  below  have  been  obtained  on 
jCDC7500  computer  with  meshes  sufficiently  fine  that 
numerical  errors  are  believed  to  be  of  minor  importance, 
is  little  scope  for  establishing  grid-independence  of 
dimensional  flows  by  successive  mesh  refinement  with  a  computer  of 
this  size  and  accordingly  our  attention  is  directed  first  at  the 
laminar  flow  results  of  Enayet  et  al  09*3)  for  R/P  ■  7.P  at  a 
Reynolds  number  of  1093.  The  computed  results  have  been  obtained 
with  20  nodes  in  the  radial  (r)  and  circumferential  (6)  directions 
and  with  100  streamwise  planes,  70  of  which  are  in  the  <»0o  bend 
itself  and  the  remainder  in  the  approach  and  downstream  tangents. 
A  preliminary  run  was  made  of  developing  flow  in  a  straight  pipe, 
the  flow  being  allowed  to  develop  until  best  agreement  was 
obtained  with  the  90°  bend  data  at  the  station  0.58  diameters 


-.i  a 
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upstream  of  the  bend.  The  computed  profiles  obtained  in  that  run 
■2 .4  diameters  upstream  of  that  position  were  then  used  as  the 
starting  profiles  in  the  calculation  proper  (i.e.  at  a  position  3 
diameters  ahead  of  the  90°  bend) . 

|  The  development  of  the  measured  and  computed  axial  velocity 

profiles  is  shown  in  figure  2  at  three  stations.  Due  to 
differences  in  the  refractive  index  of  the  perspex  tube  walls  and 
the  fluid,  the  LDA  measurements  were  made  along  the  non-parallel 
lines  shown  in  the  inset;  the  corresponding  computed  velocities 
along  the  same  lines  have  been  obtained  by  interpolation. 
Agreement  between  measurement  and  calculation  is  generally  very 
close.  This  set  of  data  has  previously  been  predicted  by  Levy  et 
al  (1983)  using  their  scheme  for  correcting,  in  a  one-dimensional 
marching  manner,  the  potential  flow  pressure  distribution. 
Agreement  between  their  results  and  the  present  computations  is 
close  over  the  first  half  of  the  bend,  but  deteriorates  somewhat 
towards  the  flow-exit,  presumably  because  their  rather  simple 
pressure  treatment  gradually  loses  accuracy  as  the  shear  flow 
becomes  more  tangled.  In  fact.  Levy  et  al  (19P3)  report  two 
regions  of  flow  separation  in  their  computations,  one  on  the 
outside  of  the  bend  from  0-15°  of  arc  and  the  other  on  the  inside 
from  about  80°  to  the  downstream  limit  of  computations  (the 
location  of  this  final  plane  was  not  reported)  .  There  seems  no 
suggestion  in  the  experimental  data  that  separation  occurred  and 
in  the  present  computations  the  minimum  value  of  friction  factor 
was  57%  of  that  in  fully  developed  flow  in  a  straight  tube  at  the 
same  Reynolds  number  and  occurred  just  beyond  the  90®  position. 

The  corresponding  turbulent  flow  calculations  are  shown  in 
figures  3-5,  the  experimental  data  again  being  due  to  Enayet  et  al 
(1982)  obtained  in  the  same  curved  duct  but  now  at  a  Reynolds 
number  of  4.3  x  10**.  Computations  of  the  velocities  were  made  on 
a  20  x  27  mesh  in  the  cross-sectional  plane,  the  additional  seven 
nodes  being  placed  in  the  'parabolic  sublayer'.  At  the  node 
adjacent  to  the  wall,  a  typical  value  of  the  normal  distance  i/UT /  v 
was  3;  thus  turbulent  stresses  were  entirely  negligible  compared 
with  those  due  to  viscous  shear.  As  with  the  laminar  flow  tests, 
a  preliminary  run  was  made  in  a  straight  duct  to  develop  the 
boundary  layer  by  the  same  amount  as  the  experiments  indicate. 
Streamwise  velocity  profiles  are  shown  in  figure  3.  For  this 
flow,  agreement  beyond  30°  is  not  as  complete  as  for  laminar  flow. 
Although  the  general  character  of  the  development  is  reproduced, 
with  the  accumulation  of  low  momentum  fluid  on  the  inside  of  the 
bend,  the  detailed  behaviour,  particularly  the  double  peak  in  the 
streamwise  velocity  on  the  plane  of  symmetry  at  68°  and  75°,  is 
not  correctly  predicted.  (In  this  region  the  present  computations 
give  nearly  the  same  behaviour  as  those  of  Levy  et  al  (1983)  who 
have  also  examined  this  flow) .  At  one  diameter  downstream  of  the 
bend,  agreement  is  significantly  improved.  indeed,  while  the 
remarks  above  have  emphasized  differences  between  experiment  and 
computation,  it  may  be  noted  that  the  level  of  agreement  achieved 
is  far  superior  to  that  obtained  by  Chang  et  al  (1983)  in  a  square 
sectioned  180°  bend  with  the  same'  turbulence  model  in  the  main 
flow  but  where  wall  functions  were  used  to  apply  near-wall 
boundary  conditions. 

Typical  secondary  velocity  profiles  are  shown  in  figure  4. 


The  'most' 'str  iking"  feature  "Is  “the- fact" that'  the  maximum  secondary 
velocity  occurs  very  close  to  the  wall,  well  within  the  region  of' 
substantial  viscous  influence.  It  is  this  feature  that  makes  the: 
wall-function  approach  inappropriate  to  this  type  of  flow.  Thet 
streamwise  wall  shear  stress  around  the  inner  line  of  symmetry  is 
shown  in  figure  5.  Due  to  the  combined  effect  of  secondary  flow; 
and  the  adverse  pressure  gradient  encountered  on  exit  from  the- 
bend,  the  shear  stress  falls  to  zero  (indeed  becomes  marginally 
negative)  in  the  exit  region.  The  present  semi-elliptic  scheme 
can  strictly  not  cope  with  streamwise  separations  but  the  reverse 
flow  velocities  were  very  weak  and  confined  within  the  viscous 
region  so  that  no  special  measures  were  needed  since  the, 
troublesome  convective  terms  were  negligible.  That  the  turbulent 
flow  should  separate  but  the  corresponding  laminar  flow  should  not 
is  of  course  entirely  contrary  to  experiences  in  two-dimensional 
flow.  j 

The  thermal  energy  equation  has  also  been  solved  in  this  case 
for  a  molecular  Prandtl  number  of  0.7.  A  uniform  heat  flux  is 
applied  starting  3  diameters  upstream  of  the  bend,  the  flow 
entering  at  uniform  temperature.  As  the  flow  develops  around  the 
bend  there  is  a  moderate  increase  in  the  average  level  of  the 
Nusselt  number.  Figure  f)  indicates  that  the  circumferential 
variation  is  small  over  most  of  the  perimeter,  but  by  PP°  a  strong 
jdecrease  in  heat  transfer  coefficient  develops  on  the  inside  of 
ithe  bend  due  to  the  fact  that  the  fluid  arriving  there  has  been 
passing  close  to  the  pipe  wall  as  it  moves  from  the  outside  to  the 
inside  of  the  bend  (cf.  figure  4).  There  is  also  a  small  peak  in 
Nu  at  the  outer  line  of  symmetry,  in  this  case  due  to  the 
jimpingement  of  relatively  cool  high  velocity  fluid.  The  ratio  of 
maximumiminimum  Nusselt  number  grows  to  nearly  5:1  by  75°.  These 
(large  variations  diminish  rather  slowly  on  entry  to  the  straight 
downstream  tangent,  the  ratio  of  outer  to  inner. heat  transfer 
coefficients  having  fallen  only  to  4:1  at  one  diameter  beyond  the 
end  of  the  bend. 

An  even  stronger  impression  of  the  acrobatic  variation  of 
heat  transfer  coefficient  is  conveyed  by  figure  7  which  shows  the 
Nusselt  number  variation  around  the  inner  line  of  symmetry.  (The 
values  are  normalized  by  the  maximum  value  in  the  field  which 
occurs  at  the  initial  station.  Since  the  entering  fluid 
temperature  is  uniform,  the  actual  value  of  Nusselt  number  is  of 
no  physical  significance  -  it  depends  simply  on  the  distance  from 
the  pipe  wall  of  the  near-wall  node)  .  The  modest  rise  in  Nu  at 
entry  to  the  bend  arises  from  the  initial  flow  acceleration  on  the 
inside  of  the  bend;  thereafter  the  strong  decrease  indicates  the 
accumulation  of  sluggish  low-momentum,  heated  fluid  on  the  inside 
transported  there  by  the  secondary  flow.  The  whole  pattern 
closely  resembles  the  variation  of  wall  shear  stress  shown  in 
figure  5,  save  that  the  Nusselt  number  does  not  go  to  zero  with^w. 
The  minimum  in  Nu,  occurring  at  the  bend  exit  is,  however,  less 
than  20%  of  that  at  entry  to  the  bend.  j 

At  present  no  measurements  of  local  heat  transfer' 
coefficients  are  known  for  developing  flow  around  pipe  bends. 
Seban  and  McLaughlin  ( ) 3 >  have,  however,  reported  such  data  for 
(the  flow  of  water  in  long  coils  where  development  lengths  are 
sufficient  for  fully-developed  flow  to  be  established.  This  case 


'has'- been  ~  simulated  using  only  four  streamwise  pressure  planes, 
downstream  values  being  successively  transferred  upstream  until 
the  flow  pattern  ceased  to  change. *  The  resultant  circumferential' 
distributions  of  Nusselt  number  are  compared  with  the 
Seban-McLaughl in  data  in  figure  8.  The  case  of  R/d  =  5?  is  quite 
well  predicted  though  the  experiments  display  a  somewhat  greater 
circumferential  variation  than  do  the  computations.. 
Qualitatively,  one  would  expect  this  type • of  difference  to  arise 
with  a  Boussinesq  viscosity  model  (BVM) .  Real  turbulence  is  known 
to  be  highly  sensitive  to  small  amounts  of  streamline  curvature, 
heat  transfer  coefficients  being  augmented  in  boundary  layers 
developing  on  concave  surfaces  and  damped  on  convex  ones.  BVMs  do 
not  capture  this  sensitivity  and  it  is  to  this  failure  that  we  may 
attribute  the  differences  shown  in  figure  Pa.  A  similar  behaviour 
is  also  shown  in  figure  Pb  for  a  bend  with  R/d  *  P.5.  The 
differences  in  Nusselt  number  between  computation  and  experiment 
on  the  outer  half  of  the  bend  are  now  more  marked,  the  computed 
levels  being  nearly  40%  below  the  data  in  places.  This  result  may 
reflect  simply  the  deterioration  in  the  model's  ability  to  mimic 
the  turbulent  transport  processes  as  curvature  effects  become 
progressively  stronger.  It  shoul.d  at  least  be  noted,  however, 
that  the  experimenters  had  several  difficulties  with  this  tighter 
bend.  The  stainless  steel  coil  was  fabricated  from  a  straight 
tube  and,  in  bending,  creasing  occurred  on  the  inside  of  the  bend, 
stretching  on  the  outside  (with  attendant  modification  of  the 
tube's  electrical  resistivity)  and  some  distortion  of  the 
cross-section  shape.  -In  view  of  these  departures  in  the 
experiment  from  the  idealized  geometry  and  boundary  conditions, 
firm  conclusions  on  this  issue  cannot  be  reached. 


4.  Conclusions 

•  A  semi-elliptic  solving  procedure  combined  with  the  parabolic 
sublayer  treatment  adjacent  to  the  pipe  wall  offers  a 
promising  route  for  the  numerical  analysis  of  complex 
three-dimensional  flow  in  pipe  bends  with  strong  curvature. 

•Although  agreement  is  not  complete,  broadly  satisfactory 
correspondence  has  been  demonstrated  with  the  velocity  data 
of  Enayet  et  al  (1982)  obtained  in  a  tight  90°  bend. 

•  The  present  Boussinesq  viscosity  model  leads  to  an 
underestimate  of  the  circumferential  variation  in  heat 
transfer  coefficient  in  fully  developed  flow  through  a  coil 
(though  for  many  engin  ering  purposes  the  measure  of 
agreement  achieved  would  be  satisfactory) .  The  differences 
are  believed  to  arise  from  the  known  insensitivity  of  BVMs  to 
streamline  curvature.  The  introduction,  in  the  future,  of  a 
turbulence  model  based  on  second-moment  closure  should  lead 
to  more  satisfactory  predictions. 


+  This  reduction  In  the  size  of  the  pressure  array  allowed  us  to  increase 
the  number  of  radial  nodes  to  40  and  the  number  of  circumferential  nodes 
to  25.  This  refinement  from  the  original  27  x  20  mesh  gave  levels  of 
HU  which  differed  by  at  most  1%  from  those  obtained  with  the  coarser  mesh, 
figure  lib. 
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Figure  1  Toroidal  coordinate  system 
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Abstract 

The  paper  points  out  that,  in  the  numerical  computation  of  elliptic 
or  three-dimensional  turbulent  flows,  the  neglect  of  pressure-variations 
across  the  very  thin  viscosity  affected  region  near  the  wall  allows  a 
fine-grid  analysis  of  this  sublayer  without  prohibitive  penalties  in 
core  or  computational  time.  The  scheme  has  been  successfully  applied 
to  the  three-dimensional  flow  around  a  U-bend. 


Introduction 


In  the  numerical  study  of  complex,  two-dimensional  turbulent 
flows  near  walls,  one  commonly  finds  that  a  different  approach 
to  handling  the  near-wall  low-Reynolds-number  region  (i.e.  the 
viscous  sublayer  and  "buffer"  layer)  is  adopted,  depending 
upon  whether  the  flow  as  a  whole  is  of  boundary-layer  or 
"recirculating"  type.  In  the  former  case,  because  an 
economical,  once-through  marching  solution  can  be  applied, 
the  near-wall  zone  is  often  analysed  by  adopting  a  fine 
grid  to  cover  the  low-Reynolds-number  region.  Such  an 
approach  is  rarely  practicable  in  an  elliptic  flow,  however, 
because  the  coupling  of  the  velocity  and  pressure  fields 
requires  an  iterative  solution;  this  feature  means  not  only 
that  computer  times  may  typically  be  two  orders  of  magnitude 
greater  than  for  a  boundary-layer  study,  but  also  that  the 
dependent  variables  must  be  stored  over  the  whole  solution 
domain.  .Because  of  the  substantial  core  and  computer  time 
requirements,  the  near-wall  region  is  usually  handled  by  way 
of  wall  functions  (1)  in  which  wall  adjacent  nodes  are 
placed  relatively  far  from  the  surface  so  that  they  lie  in 
the  region  of  fully  turbulent  fluid.  The  wall  functions 
attempt  to  embody,  through  a  mixture  of  analysis  and  experiment¬ 
al  data,  the  integrated  effects  of  the  near-wall  sublayer;  in 
this  way,  no  substantial  near-wall  mesh  refinement  is  needed. 
This  feature  is  crucial  in  keeping  the  overall  core  and 
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computing  time  requirements  to  manageable  levels  and  is  the 
reason  why,  despite  the  crudeness  of  the  physical  treatment, 
wall  functions  are  nearly  universally  adopted.  The  above 
remarks  apply  with  even  greater  force  to  three-dimensional 
flows. 

The  problem  with  such  a  simple  approach  to  the  physics  is 
that  it  is  not  adequate  to  account  for  the  diversity  of  the 
phenomena  displayed  by  turbulent  flow  near  walls.  There 
are,  for  example,  many  situations  where  the  velocity  vector 
parallel  to  the  wall  undergoes  strong  skewing  across  the 
low-Reynolds-number  region,  a  feature  which  no  wall-function 
approach  appears  to  have  mimicked  successfully. 

The  purpose  of  the  present  note  is  to  recommend  a  new 
numerical  practice  Chat  facilitates  the  use  of 
fine  near-wall  mesh  in  computing  elliptic  and  three-dimension¬ 
al  flows.  This  development  thus  opens  the  way  to  more  refined 
modelling  of  the  physics  of  the  near-wall  region  than  has 
hitherto  been  employed. 
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The  PSL  Scheme  and  its  Application 

The  PSL  scheme  is  based  on  the  idea  that,  while  the  flow  as 
a  whole  must  be  regarded  as  elliptic,  there  is  a  thin 
parabolic  sublayer  (whence  the  acronym)  immediately  adjacent 
to  the  wall  across  which  static  pressure  variations  are 
negligible  or,  in  the  case  of  highly-curved  surfaces,  where 
the  variation  may  be  obtained  by  assuming  radial  equilibrium. 
This  parabolic  sublayer  is  taken  to  extend  over  the  whole  of 
the  low-Reynolds-number  region  where  the  turbulent  transport 
properties  exhibit  such  a  strongly  non-linear  variation. 

If,  as  we  have  argued  is  desirable  on  physical  grounds,  a 
fine  grid  treatment  is  employed  across  this  region,  then 
major  simplifications  may  be  made  to  the  conventional, 
incompressible  elliptic  treatment,  (2).  Our  own  implement¬ 
ation  of  the  idea  has  been  within  the  context  of  finite 
volume  procedures  employing  a  staggered  arrangement  of 
dependent  variables,  Fig.  1  (analogous  simplifications  can 
clearly  be  adopted  with  an  orthodox  finite  difference 
method).  Within  the  PSL: 

(i)  the  pressure  does  not  require  storing  (it 
is  given  by  the  pressure  just  outside  the 
region) ; 

(ii)  thus,  no  Poisson  or  pressure-perturbation 
equation  has  to  be  solved; 


(iii)  .the  velocity  component  normal  to  the  wall  may  be 
obtained  very  rapidly  by  cell  continuity  rather 
than  by  solving  the  normal  momentum  equation. 

Thus,  referring  to  Fig.l,  for  a  Cartesian 
mesh  and  a  two-dimensional  flow, 

V(I,J)  -  (U(I-l,J)-U(I,J))6y/<5x  +  V(I,J-1) 

These  are,  of  course,  the  classic  boundary-layer  simplifications 
known  for  80  years;  what  appears  to  be  novel  is  their  application 
to  a  very  thin  sublayer  in  a  shear  flow  that  is  overall  not 
analysable  under  the  boundary-layer  approximation.  The  adaptations 
required  to  most  elliptic  solving  schemes  to  incorporate  the  PSL 
treatment  will  be  trivial.  In  the  codes  used  at  UMIST  the 
momentum  equation  (or  equations)  for  the  component (s)  parallel 
to  the  wall  are  solved  simultaneously  over  both  the  elliptic 
region  and  the  PSL;  the  velocity  component  normal  to  the  wall 
within  the  PSL  is  found  next  by  applying  continuity  to  the  pressure 
cells;  thereafter,  the  momentum  equation  for  this  component  is 
solved  over  the  elliptic  region  only.  Finally  the  pressure  or 
pressure-correctio'n  equation  is  solved  over  the  elliptic  region 
with  corresponding  adjustments  also  being  made  to  the  velocity 
field. 

There  is  also  often  scope  for  reducing  storage  associated  with 
velocity-component  information.  In  many  cases,  at  the  expense 
of  somewhat  more  code  reorganization,  the  solution  can  be 
arranged  so  that  velocities  in  the  PSL  are  stored  only  on  the 


row  contiguous  with  the  elliptic  region  and  along  two 
1-dimensional  columns  of  nodes  (which  are  successively  over¬ 
written),  rather  than  in  a  full  2-dimensional  array. 

A  referee  has  queried  the  use  of  the  PSL  approach  in  the 
vicinity  of  a  stagnation  point  where  the  variation  of  pressure 
normal  to  the  wall  is  relatively  rapid.  Perhaps  the  first 
thing  to  emphasize  is  that  any  errors  associated  with  pressure 
variations  across  the  buffer  region  will  affect  coarse-grid 
wall-function  schemes  at  least  as  much  as  a  PSL  approach. 

Our  experience  at  UMIST  suggests  in  fact  that  even  on  the 
axis  of  an  impinging  jet  the  PSL  approximation  can  be  applied 
over  most  of  the  low  Reynolds  number  region.  Although  in 
our  applications  to  date  the  PSL  treatment  has  been  applied 
to  the  same  number  of  cells  in  each  column  (viz  Fig.l),  this 
practice  is  neither  necessary  nor  optimal  in  some  applications. 
In  the  impinging  jet,  for  example,  a  thin  PSL  at  the  stagnation 
point  could  be  expanded  to  cover  the  full  height  of  the  domain 
once  the  jet  had  been  deflected  into  a  radial  wall  jet.  Indeed 
a  self-adjusting  scheme  for  the  number  of  nodes  in  the  PSL  at 
any  x-position  could  readily  be  devised. 

The  most  important  field  of  application  of  the  approach  is 
perhaps  in  three-dimensional  flows  describable  by  the  partially 
parabolic  equations  for  there  only  the  pressure  field  requires 
3D  storage.  The  scheme  has  been  successfully  applied  by  the 
authors  to  the  turbulent  flow  in  a  circular  tube  around  a  90° 


each  radial  strl 

bend.  Mine  nodes  have  been  put  in  the  parabolic  sublayer  <along/  nod 
For  identical  grid  densities  in  the  fully  turbulent  region, 
computing  times  are  no  longer  (in  fact  somewhat  less)  than 
with  our  previously  used  wall-function  approach;  core  require¬ 
ment  is  also  little  affected  because  most  of  this  is 
associated  with  the  pressure  field  which  (alone)  has  to  be 
held  on  a  three-dimensional  array  (and  which  is  identical 
for  the  two  approaches  since  there  are  no  pressure  nodes 
within  the  PSL).  The  PSL  scheme  has  also  been  applied  in 
high  Reynolds  number  lamina T  flow  around  pipe  bends  where 
again  near  the  wall  the  velocity  field  undergoes  very 
rapid  changes  but  the  pressure  is  obtained  adequately  via 
radial  equilibrium. 

The  scheme  has  since  been  adopted  by  two  of  our  colleagues 
who  had  hitherto  been  using  a  fine-grid  low-Reynolds-number 
approach  within  two-dimensional  fully  elliptic  treatments. 

When  there  were  no  flow  reversals  in  the  near-wall  layer, 
the  introduction  of  PSL  reduced  their  computational  times, 
in  one  case  by  a  factor  of  two  and  in  the  other  by  a  factor 
of  three ,-t  Benefits  were  much  reduced  when  reverse  flow 
was  present  but  it  seems  likely  that  these  can  be  substantially 
restored  by  reorganizing  the  solution  in  the  PSL  so  that  the 
direction  of  marching  is  always  that  indicated  by  the  velocity 
at  the  outer  edge  of  this  sub-layer. 

t  Their  numerical  results  obtained  with  PSL  are  insignificantly 
different  from  those  given  by  the  fully  elliptic  solution 
using  an  identical  mesh. 


Conclusion 


The  PSL  approach  allows  a  fine-grid  resolution  to  be  applied 
to  the  near-wall  sublayer  with  no  significant  increase  in 
computer  time  or  storage,  compared  with  a  conventional  wall- 
function  treatment.  Because  the  former  scheme  facilitates 
a  better  modelling  of  the  turbulent  transport  processes,  it 
is  thought  that  in  many  cases  it  may  supplant  the  latter. 
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